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SUMMARY 


An analysis has been developed and a computer code, P465 Version A, written for 
predicting three-dimensional transonic potential flow about inlets, ducts, and bodies. 
The basic approach is very general and essentially independent of geometry. The 
analysis, as programmed, uses cylindrical coordinates and admits boundary 
conditions for the aforementioned configurations. Solutions of the full potential 
equation are obtained using finite differences and successive line over-relaxation 
(SLOR). Extrapolation and a sequence of increasingly denser grids are used to 
accelerate convergence and thus decrease the computer time required to obtain a 
solution. 

The analysis has been programmed in extended FORTRAN IV for the Control Data 
Corporation CYBER 203 computer. Use of reasonably dense meshes for 
three-dimensional computations requires either a very large computer, such as the 
CYBER 203, or a great amount of complex software to use disk memory for backup of 
core memory. The run time on the CYBER 203 for a t3rpical inlet calculation is about 
five minutes. Some special features of the CYBER 203 have been used in the 
computer code to maximize efficiency of the code. Use of the analysis on another 
super computer, such as the CRAY 1, is practical but would require a certain amount 
of recoding. ' 

Comparisons between computer results and experimental measurements are 
presented to verify the capabilities of the code. Agreement between analysis and 
experiment is excellent for flowfields which are essentially inviscid and irrotational. 
The greatest disagreements between the analysis and experiment occur for flows 
where viscous (boundary-layer) effects are significant. 

Program input and output formats are described and examples are presented. 

This report does not discuss geometry specification. The analysis code requires input 
of a computational mesh and all intersections of the computational mesh with 
surfaces, including coordinates and components of the surface normal. This geometry 
specification can come from any source and the file of intersections can be in any 
order. The program does extensive checking of the geometrical input. 
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1.0 INTRODUCTION 


This report describes a computer code, P465 Version A, for predicting the inviscid 
irrotational flowfield about a three-rdimensional inlet, duct^or body. The flowfield is 
generated by a finite-difference, line-relaxation solution of the complete potential 
flow equation. The predicted results are excellent approximations of actual flowflelds 
where viscous effects are small, i.e., flows with thin attached boundary layers. 

During the past decade there has been tremendous progress in the field of transonic 
potential-flow computation. The work started with Murman and Cole^^^ who used 
upwind differencing in supersonic flow regions to solve mixed elliptic and hyperbolic 
equations. Their work was based on small disturbance theory where the flow was 
principally aligned with one coordinate. Jameson^^^ extended this approach to the full 
potential equation with arbitrary flow direction. Several classes of three-dimensional 
geometries have been solved using the full potential equation^2-5) Most extensions to 
three dimensions have been based on some mapping from the physical coordinates to 
a more convenient computational coordinate system. If such mappings can be found, 
the computational problem, exclusive of the generation of the mapping, can be 
greatly simplified. This approach becomes more difficult as the geometry of objects to 
be considered becomes more complex, because it becomes more difficult to find 
convenient mappings that will provide reasonable coordinates for the computation. A 
further disadvantage with a complex mapping is that it can be very difficult to have 
any physical understanding of the solution process. 

The approach used for this analysis separates the coordinate system and the 
geometry. Thus, this approach does not require the finding of a new transformation 
for every family of geometries. The method was originally tested for axisymmetric 
problems^®^ and later applied to axisymmetric inlets at angle of attack^^\ This report 
describes the extension of the method to general three-dimensional geometries. The 
method also has the advantage that the solution is obtained in a physical coordinate 
system with physical variables and hence the entire solution process is more 
understandable. The analysis, as programmed, uses cylindrical coordinates since they 
are suited to nacelle calculations, but the analysis could easily be changed to use a 
Cartesian, skewed Cartesian, or any other coordinate system. 

The analysis is independent of the coordinate system except that if the coordinate 
system is a poor choice for the problem, accuracy will suffer as there will be 
insufficient mesh nodes. The practical limit on the number of mesh nodes is dictated 
by computational time. The nature of the relaxation process imposes an exponential 
increase in computational time with increasing mesh density. This severely limits the 
number of mesh nodes and makes it essential that the coordinate system be such that 
the mesh nodes are well placed. The cylindrical coordinate system used in the 
computer code gives reasonable node placement for nacelle, duct, and body geometries 
where such geometries are approximately bodies of revolution. Cylindrical 
coordinates can be used for geometries that are not nearly bodies of revolution, but 
they may not be very efficient. 



2.0 ANALYSIS 


Solutions to the full partial differential equation for compressible potential flow along 
with appropriate boundary conditions are desired for three-dimensional flow. 
Approximate solutions of this equation are obtained by replacing the partial 
differential equation and boundary conditions with finite-difference equations for the 
value of the potential function, at discrete points in the flowfield and on the 
boundary. The finite-difference equations are obtained by replacing partial 
derivatives in the differential equations with difference quotients. The difference 
equations are coupled, nonlinear, algebraic equations in the value of the potential 
function at field and surface mesh points. The nonlinear difference equations are 
solved by an iterative process to the desired accuracy. The degree of approximation of 
the difference solution to the true solution of the partial differential equation is a 
function of the mesh spacing. 

The solution technique can use an arbitrary mesh as long as the mesh is dense 
enough to ensure that there are several mesh nodes between independent surfaces. As 
programmed, a cylindrical mesh with variable mesh-line spacing is used because it is 
relatively efficient for inlet and nacelle configurations. A typical constant 6 cut 
showing the mesh and an inlet cowl is displayed in figure 1. The mesh is not required 
to intersect with the body in any predescribed manner. Nodes used in the calculation 
are mesh intersects in the flowfield and intersections of mesh lines with surfaces. 
Difference quotients are written to take into account variable spacing between field 
and surface nodes. 

The difference equations consist of a large number of coupled, nonlinear, algebraic 
equations. The direct solution of such a set of equations is certainly impractical if not 
impossible. A multistep process is used to develop a solution procedure for the 
difference equations. The basic approach uses an iterative technique, starting with an 
initial guess, computing an approximate solution, and then repeating the cycle or 
sweep with the approximate solution used to compute a better approximate solution. 
As a minimum requirement, the iteration technique has to be convergent. In practice, 
for the three-dimensional transonic potential-flow problem, the iteration technique 
has to converge rapidly enough that the computation time is practical. 

The first step in the iterative solution process is to linearize the difference equations. 
This is done by dividing terms into coefficients and unknowns, and using the initial 
guess or solution for the previous iteration to approximate the coefficients. The linear 
equation system that results is fully coupled, or implicit, and would have to be solved 
simultaneously. Simultaneous solution of such a large system of linear equations is 
not practical. As a consequence, the equations have to be at least partly decoupled by 
substituting old values of the potential function, </>, for some of the (unknown) values 
of (j) (denoted by </>+). How this substitution is done affects stability, convergence, 
convergence rate, and the practicality of the solution procedure. The process of mixing 
old and new values and iterating (sweeping) to obtain a solution is commonly called 
relaxation and is discussed in many references. (See ref. 7, for example.) Generally, 
the more thoroughly the equations are decoupled, the less work is required to take 
one relaxation sweep,* however, the more relaxation sweeps are required to reach a 
given level of convergence or degree of accuracy of the solution. A fully decoupled 
system, where only one equation needs to be solved at a time, is known as explicit 
and generally exhibits very slow convergence properties. 
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Figure 1. Typical Coarse Computational Mesh for a Constant d Cut in the Vicinity 
of the Inlet 
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A decoupling of the equation system so that all the equations along one mesh line are 
solved simultaneously is known as line relaxation. It is a popular technique because 
there exists a very efficient algorithm known as the Thomas or tridiagonal algorithm 
to obtain the solution for such an equation system. The relaxation process can be 
accelerated or stabilized, as desired or required, by processes known as over- or 
under-relaxation, respectively. Over-relaxation consists of taking the difference of the 
new and previous solutions, multiplying the difference by a factor between one and 
two, and adding the difference to the previous solution to create an improved new 
solution. Under-relaxation is similar except the multiplication factor is less than one. 
In this code over-relaxation is used for one term of the partial differential equation 
and under-relaxation for another. The combination of line relaxation and 
over-relaxation is commonly called successive line over-relaxation (SLOR). This is 
basically the technique used in the analysis discussed herein. 

A penalty involved in using a line-relaxation scheme on a vector computer is that the 
tridiagonal algorithm, used to solve the simultaneous equations for a line, does not 
vectorize efficiently. However, the ratio of computing time spent solving the 
equations with the tridiagonal solver, versus the time spent calculating the equation 
coefficients, is small. Also, the SLOR scheme takes many less relaxation sweeps to 
reach a given level of convergence than an explicit technique. A sweep is a complete 
solution in sequential order of all the difference equations for the flowfield. 

Each relaxation sweep is computed using SLOR. Additional techniques have also 
been employed to accelerate convergence. One such technique is extrapolation using 
intermediate solutions. Extrapolation is based on the theory of linear algebraic 
equations and the existence of eigenvalues and eigenvectors. When the solution 
process appears to be developing in a particular manner it can be extrapolated to 
obtain a much better approximation of the final solution. Another technique is the 
use of a sequence of meshes of increasing density. The initial field is assumed to be 
uniform flow and is applied on a very coarse mesh. Relaxation is continued on this 
mesh until partial convergence is obtained. The initial mesh is too coarse to be used 
to compute an accurate solution to the flow problem, but the solution from the coarse 
mesh is a much better starting guess for the next mesh than uniform flow. The 
advantage of using a sequence of meshes is that computation is cheaper per sweep, 
and convergence is faster on a coarser mesh. 

Along with the consideration of the practicality of making a relaxation sweep, the 
linearization and decoupling of the difference equations has to be such that the 
relaxation process converges. One aspect critical to the transonic flow solution is the 
use of upwind differencing where the flow is supersonic. 

The organization of this analysis, as programmed, is biased by the strengths and 
weaknesses of the computer for which the code was originally written. In particular, 
the CDC CYBER 200 family of computers has large memory and operates most 
efficiently with long- vector operations. The large memories make three-dimensional 
codes possible. The efficiency of long-vector operations makes it desirable to decouple 
the equations so that many coefficients can be calculated at once. 



2.1 EQUATION AND BOUNDARY CONDITIONS 


The equation to be solved is the complete equation for inviscid, irrotational flow 
formulated in terms of a velocity potential, <j>. This form of the equation has the 
significant advantage that there is only one unknown that has to be stored at mesh 
nodes. This is particularly important for a three-dimensional elliptic problem where 
computer storage is very critical. 

2.1.1 POTENTIAL FLOW EQUATION 

The equation for the velocity potential in cylindrical coordinates is 




( 1 ) 


where </> is the velocity potential and the local speed of sound, a, is given by 

2 2 Y-1 / V 2 ^ / 2 ^ 2\ 

a = " 2 

The potential equation in Cartesian coordinates is 


(2) 
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This form of the equation is used in generating difference equations for points on the 
axis. 


The velocity components in the flowfield (u, Uj., u^) for the cylindrical coordinate 
system are obtained from the potential function with the following relationships: 
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The velocity components for a Cartesian reference frame (u, v, w) are obtained using 
the equations, 


V = ^ = u cos 0 - sin 0 

y r 0 

w = / = u sin 0 + u_ cos 0 . (5) 

z r 0 


2.1.2 BOUNDARY CONDITIONS 

The boundary condition at solid surfaces is that the velocity normal to the surface, 
equals zero. The boundary condition at the exit of a duct, or at the compressor face of an 
inlet, is that the axial velocity is fixed at the uniform value that gives a specified mass 
flow. At the left boundary of the computational field, the potential function, <f>, is 
specified. For an inlet flowfield computation (figure 2), at the left boundeiry. 


/ = u^x + v^r cos 0 + sin 0 + constant. 


( 6 ) 


This is also the equation for the initial flowfield approximation. 

For an isolated inlet in a freestream, the following additional boundary conditions 
are used. Equation 6 is used to specify <f> for 90 degrees ^ 0 ^ 270 degrees and 
r 'I'he outflow velocity, 4>^, for the computational cylinder is specified by 




r=r 


max. 




r=r 

max. 




(7) 


for 0 degrees ^ < 90 degrees and 270 degrees < 6 ^ 360 degrees and r = . 

is specified as equal to u^ at the right boundary of the computational field outside 
of the inlet. 

The boundary conditions on the far field for an isolated nacelle or body are only 
approximate. The correct boundary condition is that the velocities approach the 
freestream values far from the body. This boundary condition cannot be directly satisfied 
for a potential-flow solution in a finite computational volume. As written, the boundary 
conditions specify either the normal or tangential component of velocity at any surface. 
In most cases the normal component of velocity is specified for outflow surfaces and the 
tangential component for inflow surfaces. The exceptions are nonzero yaw angle or 
negative angle of attack. For these cases, the boundary conditions on some of the 
surfaces could be changed; for example, the analysis is convergent for a negative angle of 
attack which essentially reverses the boundary conditions on the top and bottom of the 
computational cylinder. Changes to the far-field boundary conditions would probably 
affect accuracy and convergence rate of the solution, however, no extensive studies of 
these effects have been made. 




Figure 2. Geometry and Boundary Conditions for iniet Fiowfieid Computation 



The present far-field boundary conditions are satisfactory provided the boundaries are 
located far enough from the body. Locating the far-field boundaries five to ten 
nacelle diameters away from the body is usually sufficient to ensure negligible effects of 
the boimdary conditions on the computed flowfield near the body. If the far-field 
boundary is so placed, the effects of the boundary conditions are of the order of or 
smaller than other approximations such as treating a nacelle as an isolated nacelle. 

2.2 WEIGHTING AND ALTERNATE FORMULAS 

A possible problem with this approach to solving the potential flow equation is that the 
step size between a surface point and an adjacent field point is uncontrolled and can 
become arbitrarily small. Derivatives are calculated by dividing by the step size and thus 
any error in the potential function, <f), can be magnified by the reciprocal of the step size, 
which can become arbitrarily large. This can cause accuracy and stability problems if not 
considered. To avoid this problem, alternate formulas that do not use the field point </> 
value are used when the step size to a surface point is very small. To obtain a smooth 
solution process, a weighted combination of the regular and alternate formulas is used if 
the surface points are closer to the field points than one-half the local field-point spacing. 
This procedure, which is called weighting in this document, has been used for many 
difference quotients. Typical formulas are presented in this report. 

The <f) values at field points adjacent to surface points are used in differencing at other 
field nodes, even when they are very close to surface nodes. It is essential that as a field 
and surface node become close, the values at these nodes approach the same value. 
This behavior is ensured by calculating the <f> value at the field node by interpolation 
between the surface value and other field-point values, if the spacing to the surface is 
very small. A weighted combination of interpolation and the usual field-point difference 
equation is used when the spacing between the field node and any surface node is less 
than one-half the local mesh spacing. 

2.3 DIFFERENCE QUOTIENTS 

Difference equations are generated by replacing partial derivatives in the partial 
differential equation with difference quotients. The choice of difference quotients is not 
arbitrary, nor is it unique. The difference quotients chosen are based on several, often 
conflicting, criteria. First, the difference quotients chosen must lead to a stable 
convergent solution process. Second, the use of more nodes in a difference quotient can 
give greater accuracy, but normally at the price of greatly increased complexity. For this 
analysis, where the body can intersect the mesh arbitrarily, the use of complex, 
many-point difference quotients would be extremely difficult. 

The highest-order derivatives in the partial differential equation are second derivatives, 
which require at least three points in the difference quotients. As any more than three 
points in the difference quotients would greatly complicate many aspects of the problem, 
three-point difference quotients have been used. 

There are three general categories of difference quotients: central, forward, and 
backward differences. Categorization is based on where the additional nodes used in the 
differencing lie with respect to the node at which the derivative is to be evaluated. 
Forward or backward differences relate to the use of additional nodes which lie ahead of 
or behind, in time or space, the node where the derivative is to be evaluated. In general, 
central differencing is preferred because of smaller truncation errors, however, in some 



situations, central-difference quotients may result in unstable difference equations. In 
particular, for this problem, stability requires the use of one-sided differences for second 
derivatives where the flow is supersonic. These differences are termed "upwind” as the 
direction of the differencing is opposite the flow direction. 

2^.1 FORMULAS FOR FIRST DERIVATIVES (VELOCITIES) 

The difference quotients for first-degree partial derivatives are central differenced for 
greater accuracy. Three-point differences are used giving second-order accuracy with 
truncation errors proportional to the square and/or product of the two step sizes involved. 
If surface points are used, alternate formulas may be weighted into the calculation to 
avoid problems due to dividing by the step size to the surface when the step size becomes 
very small. The first derivatives represent the flow velocities except for which equals 
ru^. Difference quotients and weighting formulas are given here for x derivatives. The r 
and S derivatives are calculated the same way. Topical configurations for differencing 
axe shown in figure 3. 

All old 0 values are used in the first-derivative calculations, except for the calculation of 
<f>^ in the last term of equation 1. Old <f> values are used for the first derivatives because 
they are used to calculate the coefficients of the terms in the difference equations before 
the difference equations are solved. 

For a regular field point (or a field point that is not adjacent to a surface point in the x 
coordinate) (fig. 3a) 
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For a field point that is surface-adjacent in x (fig. 3b) 
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a.) REGULAR FIELD POINT 



b.) SURFACE ADJACENT FIELD POINT 



c.) SURFACE ADJACENT IN OPPOSITE DIRECTIONS 
FIELD POINT 

Figure 3. Nomenclature for Typical Difference Quotients 
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where 


AXi /|^0 

< 1^0 

1.0 

IV 

T 

O 


</>x|g is the value of for the surface point S that is calculated from the surface-point <{> 
values according to the procedure described in section 2.4. 

If the field point is surface-adjacent in x in the opposite direction as well (fig. 3c) 


and 




( 3 ) 






p' = 




1.0 




P" = minimuia(P ,P ' ) 


9^J = 


( 12 ) 


(13) 


2.3.2 DIFFERENCE QUOTIENTS FOR </)^^, AND cl>^^ 

The difference quotients used for replacing the second derivatives affect the stability of 
the relaxation scheme. Central differencing is preferred for accuracy, but leads to 
unstable behavior if the flow becomes supersonic. The technique for avoiding this 
problem is to switch to upwind difference quotients when the flow becomes supersonic. 

It would be desirable if all new 0 values could be used in the difference quotients, but 
that leads to a coupled equation set that is impractical to solve. The equations presented 
here are coupled only along radial lines and are of tridiagonal form. There is a relatively 
simple, well-known, efficient, algorithm for performing a direct solution of tridiagonal 
systems. The algorithm is presented in appendix A. For x derivatives, new values are 
used at i-1 because they are available and an over-relaxation parameter, is used to 
speed convergence for the subsonic central-difference quotient. The central-difference 
quotient for uses old values at k-1 and k+1 in order that all coefficients for an x 
equals constant plane can be calculated at one time. This is desirable to exploit the 
vector speed of the CYBER 203. An under-relaxation parameter, o)^, is used to maintain 
stability. 
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Central-difference formulas: 



U)/^ 

0 


Z u 


naximum ^ u ) 



( 18 ) 


cojj and are reference over- and under-relaxation parameters. They are the values used 
for a uniform mesh. The values and cog have been modified for the local mesh spacing. 
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The same formulas are used if a surface point, or points, are involved in the differencing, 
except that the appropriate surface point <f> value(s) and step size(s) are substituted. 

The upwind difference formulas for <^>j^,and use a mixture of old and new <f} 
values for stability and improved rate of convergence. Where a new value of <f> is shown 
in a formula, an old value is used, if the flow direction is such that a new value is not 
available. 

The switch between central and upwind differencing takes place as the coefficient of the 
second derivative changes sign. ITiis leads to a smooth stable behavior of the solution 
process. As an example, central differencing is used for if upwind 

differencing if 1 I 

Upwind difference formulas are shown only for 0^x for one flow direction. Difference 
quotients for r and 0 and other flow directions are similar. For a point that is regular 
in x and for l^xl 




_ ~ ^1-1, j,k ) ~ ~ ^1-2, j,k) 


(19) 


XX 
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If the point (i-l,j, k) is surface adjacent in the x coordinate with the surface between 
(i-2, j, k) and (i-1, j, k), the following formulas are used: 
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I ^1^2K-*^2) 


( 20 ) 
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( 21 ) 
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where 


“ 2 / 


<(^0 < 




P = 


1.0 


^ 2>^1 


and S is the surface point and Ax 2 is the step size between the surface point and 
(i-l,j,k). 
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If the point (i, j, k) is surface-adjacent in negative x, 



1^1 


S 


( 23 ) 


where Axj^ is the step size between point S and (i, j, k). 


2.3.3 DIFFERENCE QUOTIENTS FOB <f>r0, AND 

For subsonic flow, central-difference formulas are used for accuracy. The formulas are 
obtained by setting 


4e 

^0x 


dr 

d0 

50 


(24) 


and using formulas 8 or 9 with appropriate changes. As an example, for a field point that 
is adjacent to a surface point S and lies below S, 

2 ^x 


^1^2 (^1-^2) 







B 


(OK\ 


Weighting is not necessary in these terms as the resulting difference equations are 
weighted (see Section 2.5.2). 

The difference quotients for the cross-derivatives are also switched between central and 
upwind formulas depending on flow Mach number. The switching procedure is somewhat 
more complicated because the coefficients of these derivatives do not necessarily switch 
signs an 5 Mvhere. 
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II III 


II 



The upwind-difference quotients for cross-derivatives are relatively simple, but there are 
many variations depending on flow direction. Representative examples of upwind 
differencing for a field point that is regular with and <f>p>0 (fig. 4a) are 




4 ^ - 4 ^ - 4 ^ + 4 ^ 


xr 




(26) 


For the configuration of figure 4b, 




I 


if j-lfk 



esc 



(27) 


For the configuration of figure 4c, 







CSC 


2 


Ax 


(28) 


Although, the configuration shown in figure 4d is unlikely since supersonic flow is shown 
coming out of the surface, it can happen under unusual circumstances, particularly for 
very coarse meshes. For the configuration of figure 4d, 




Ax 


Ar 



(29) 


There is a mixture of new and old values in this formula. The philosophy is to use new 
values for (f> where they are available, consistent with the line relaxation, except in the 
case of equations 27 and 29 which are of the form - </>x)/Ar. 

The code uses central differencing if the local Mach number is less than one and upwind 
differencing if either velocity component is supersonic. The upwind- and 
central-difference formulas for cross-derivatives are weighted together between these two 
conditions. In summary, if M<1, 

^xr ^xr I central differenced 
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a.) REGULAR DIFFERENCING b.) SURFACE AFFECTED 



c.) SURFACE AFFECTED d.) SPECIAL SURFACE AFFECTED 

Figure 4. Example Configurations for Upwind Differencing of Cross- Derivatives 
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and if | 03 ^[>a or krl>a 


<#»xr = <#*xr| 


upwind difference, 


otherwise, 


4 = H 

'^xr ^xr 


central diff. 


+ (1-3V, 


xr 


upwind diff. 


(30) 


where 


xr 


/ 2 2 \/ 2 2 \ 

Ca -u )(a -u^ ) 

/2 2\/2 2 \ 

(q -u )(q -u^ ) 


(31) 


It can be seen that = 1.0, when M = 1.0 and /S^r = 0 when |u| or |uj.| = a. 

2.3.4 SPECIAL (9-DIFFERENCE QUOTIENTS 


Special difference quotients for <f>g and are derived in reference 3. These formulas are 
equivalent to the usued formulas to the order of Ad^. They are considerably more 
accurate than the standard formulas if A6 is large (AO = tt/ 4 or tt/2) and <f> is of the form 


j(((x,r,0) = jz^^(x.r) + jz^gCx.r) sin 9 + cos 9 


(32) 


<l> has approximately the above form in the far field for inlet computations, and in the 
near field if the inlet is nearly axisymmetric and centered on the axis. The advantage of 
the special formulas is improved accuracy for very coarse meshes. Such very coarse 
meshes can be desirable if a sequence of computational meshes is used for convergence 
acceleration. The formulas are 


(33) 

- - °°^u) 
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(34) 




00 


^ l^sinz^e^ ^1 - cos^^^ + sinA0^ ^1 - cos^^'jlj 




The use of the above difference quotients instead of formulas 8, 9, and 17 is available as a 
program option. 


2.4 SURFACE DERIVATIVES 

Before the field is swept, the velocity is calculated at all surface points using all old <f) 
values. The velocity components and are computed first. These are the 

components of velocity along constant x, r, and 0 cuts of surfaces, and, except at end 
points of the cuts, are computed using central-difference formulas similar to equation 8. 
One-sided differences are used at end points. Weighting is used to eliminate possible 
problems due to very closely spaced points. For any given surface point, only two of these 
velocity components can be computed in this manner. As an example, 0g^ cannot be 
computed directly for an x-intersection surface point, that is, a point created by a mesh 
line parallel to the axis intersecting the surface. 

At any point the remaining surface-derivative can be computed from 






^ "0 




0 (35) 


if the magnitude of the direction cosine (n^ with (f>^, etc.) is not very small. If the 
magnitude of the direction cosine associated with the velocity component to be computed 
is very small, the above formula can be very badly behaved since it involves the 
reciprocal of the direction cosine. In such a situation, the missing velocity component is 
calculated by interpolating along one of the two cuts through the point. The cut selected 
is the one for which the interpolation will be the most accurate. 

After all three surface velocities have been computed for each surface point, the velocity 
components (^)J.,and for surface points are computed from 
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where ^ 

These velocity components are used in the calculation of the cross-derivatives, 
</)j. 0 ,and for subsonic flow, in the surface-point boundary-condition (eq. 47) under 
certain circumstances, and for the surface properties printed for program output. 

2.5 FIELD-POINT DIFFERENCE EQUATIONS 

A difference equation is generated for each field point. The difference equation is 
generated by dividing the partial differential equation (eq. 1 ) into terms and then 
splitting each term into a coefficient and an unknown. The unknown is the second 
derivative in the term, except for the last term of the equation where the unknown is 
taken as 0 ^.. The coefficients involve only first derivatives and are calculated using 
all old values of cf>. The remaining part of the difference equation is linear in <f> and 
theoretically could be solved simultaneously for the new potential values, but as 
stated earlier, this is not practical. To make the solution process simpler for a single 
sweep or iteration, the unknowns are generated using difference quotients that have 
a mixture of old and new (f) values. 

The terms in the difference equation are listed in Table 1 , including the coefficient, the 
unknown term, and an index to the difference quotients used in differencing for the 
unknown terms. 

The difference equations are generated using a nonconservative approach. The terms 
conservative or nonconservative refer to whether or not the differencing scheme 
explicitly conserves mass. The principal flow region where this nonconservative scheme 
fails to conserve mass very well is at shocks. Other than at shocks, mass conservation 
can be improved by using a denser mesh. The errors caused by a failure to conserve mass 
exactly, as determined by comparison with experimental results, appear to be relatively 
small. 

2.5.1 AXIS 

Points on the axis are a special case because the potential equation expressed in 
cylindrical coordinates is singular on the axis. This is a problem with the coordinate 
system and not a physical problem with the flow. This problem is resolved by requiring 0 
degrees, 90 degrees, 180 degrees, and 270 degrees to be in the 6 mesh, and using the 
Cartesian form of the potential equation (eq. 3) at nodes on the axis. 
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For nodes on the axis, the terms in the potential equation are given in Table 2. 
Coordinates and notation are shown in figure 5. The points A, B, C,and D can be field or 
surface nodes. The differencing for <f>j^ and <f>j^ is the same as for the equation in 
cylindrical coordinates. The differencing for 0y and <}>z equation 8 or 9 with 

appropriate </> values and step sizes. The differencing for <f)yy and <f)^ uses difference 
quotients (eq. 16) with appropriate <f> values and step sizes. 

Table 1. Terms and Difference Quotients for Field- Point Difference Equations 


TERM 

COEFFICIENT 

UNKNOWN 

FORMULA(S) FOR 
SUBSONIC FLOW 

FORMULA(S) FOR 
TRANSONIC FLOW 

REGULAR 

SURFACE 

ADJACENT 

REGULAR 

SURFACE 

ADJACENT 




1. 

2^2 

a 

^xx 

14 

14 

19 

20,21, 

22,23 

2. 

2/2 

a 

^rr 

16^ 

16^ 

19 

20,21, 

22,23 

3. 

a 

\ r ' r 

CD 

17,34** 

17,34** 

19 

20,21, 

22,23 

4. 

- 

^xr 

24* 

24* 

26 

27,28,29 

5. 

- 

r 


24* 

24* 

26 

27,28,29 

6. 

- 

X 

<r> 

24* 

i 

24* 

26 

27,28,29 

7. 

(4)r 

4 

8^ 


8^ 

9"^ 

1 - 


* USES ALL OLD VALUES OF 0 
^ USES ALL NEW VALUES OF 0 ,( 0 

** OPTIONAL FORMULA 
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where u^j^, ue|g,and «^>z|j ^ ^ have been previously calculated. 

A typical supersonic-flow difference quotient for <f>y>a and <l>j>a is 






+ 


B 


( 39 ) 


Changes are made in the differencing at the axis if flow symmetry planes are used. If the 
flow is symmetrical about the plane defined by 0 = 0 degrees or 180 degrees, equals 
zero on the plane of s 5 Tnmetry, and the potential equation for 0 = 0 degrees or 180 
degrees becomes 

(a^-^ + a^^ - Zi 4 ^ - 0 • (40) 

\ /'^xx V /^yy ^zz '^x'^y'^xy 


Y 


0 = 0 



I 

I 


Figure 5. Notation for the Axis 
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The difference quotient is obtained by setting = (f>Q and Azg = Az^^* plane 

defined hy 0 = 6^ is also a plane of symmetry, then <f)y and <j>^ must be equal to 0.0 at 
the axis, and at the axis the potential equation is 

\ /^xx ^yy 


and it can be shown that 

^ ) / 

^zz V 

where ^ = r for 0 


2.5.2 MODIFIED DIFFERENCE EQUATION FOR FIELD POINTS ADJACENT 
TO SURFACE 

The finite-difference equation for a field point next to one or more surface points can 
cause numerical problems if the spacing to adjacent surface points is small relative to 
the local mesh spacing. The analysis uses the flow equation at the field point when the 
field point is not close to surface points. An interpolation equation is used to obtain the 
field-point </> value when the surface and the field point are very close. This is 
accomplished by the following formulas which weight between the flow equation and 
interpolation formulas depending on the step size(s) between the field point and the 
surface point(s). 

The weights IBj.,and for the x, r, and 0 interpolation formulas, respectively, are 
computed by using the following relationships. The equations presented are based on 
the configuration of figure 3c. Cubic interpolation is used. 


= 0 


sin^e, 


(41) 


(42) 


Plx 


1 - 

1 0 
0 


< AXq 

if the field point is not 
adjacent in x to a 
surface point 


(43) 


The formulas for and are similar. 


^x " ^Ix " 2^1x(^lr‘^ie)‘*‘ 3^1x^lr^ie ' 


23 



The formulas for and /3^ are similar. 

The interpolation equation for x interpolation is 


K = Ti2(3_2n)^g‘^ + (i-n)^(i+2n)j!^3* 


where 


n 


- Ti(i-n) 







and 

= (^1 + ^*2) 


(45) 


and <^B* is either or <^ 3 '*' if a new value of <f> is available at B. 

The final difference equation for (f> at the field point A is 

(1 ^ ^ ^ (Difference representation of the potential equation) 

+ (x-interpolation equation) 

+Cj 33 c (r-interpolation equation) 

+ (6 -interpolation equation) = 0 (46) 

where C is the coefficient of in the regular-difference flow equation for point A. 

2.6 SURFACE-POINT DIFFERENCE EQUATIONS 

The boundary condition at solid surfaces is that there is no flow through the surface, 
which requires that the velocity normal to the surface, be zero. The velocity, 
can be expressed in terms of its components and the direction cosines (n^^, nj., n^) of 
the unit surface normal, n, 

^ = n + n V + n.— 

r r 0 (47) 

r 

Referring to figure 6 for notation, the difference equation for the surface boundary 
condition for a surface point created by a radial mesh line intersecting the surface 
(r-intersect surface point) is 



n 




+ ( 1 - 3 ) 




L “2 


n 




+ n ^ 


pii 

r 


-il 







n„ 


^ + 




2_s :z: . ^ 

A0 


p’ 



0 , 


where 



(a“i/ ^0)^ 

1 


2 ^^ < 


2^^ 


> ^ 


0 * 


( 48 ) 


r 



Figure 6. Nomenclature for Typical Surface-Point Boundary-Condition Equation 
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Currently (f>j^ is taken as zero, however, the difference equation is expressed in terms 
of anticipating a later version of the code accounting for boundary-layer 
displacement effects by using appropriate values. 


The values of (f> at the points P' and P" are obtained by cubic interpolation using the 
values of (f> and <f>y, at points Pq and P 2 . As an example 

4. = n^(3-2n)4 . + (i-n)^(i+2n)4 , 

0 ^2 

- n(i-n) 






P * 


(49) 


where 

n = ^0 


and 




= - 


The values of at P' and at P" are obtained by linear interpolation between 
points Pq and P2. As an example, 




p. 




(50) 


The point Pq may be a field or surface node depending on geometry. 


If the point P' or P" does not lie in the flowfield (i.e., is interior to the surface), 
equation 48 is modified by using the value for or <f>j^ (for P' or P" respectively) 
computed by differencing along the surface. (See Section 2.4.) 

The formulas for 0 or x surface intersections are similar. 

If the surface point is also a field-point mesh node, P' and P" are field points, point B 
is not used, and /3 = 1.0. If the surface point is on the axis or adjacent to the axis, 
some modifications to the formulas are made. 

2.7 SOLUTION PROCEDURE 

For problems such as this one which require significant quantities of computer 
resources (central-processor time, memory, etc.), it is essential that the code be 
organized for maximum efficiency. This includes such aspects as organizing the 
computation to minimize the amount of information transferred and the number of 
times information is transferred between core and disk memory. This can be done by 
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careful planning of the sequence of computations as well as reducing memory 
requirements to a minimum. Computer central-processor time requirements can be 
greatly affected by the sequence in which order independent quantities are 
calculated. 

2.7.1 CODE ORGANIZATION 

One of the major problems with three-dimensional codes is storage in the computer. 
Storing a three-dimensional field when using a tjrpical mesh requires a large number 
of words. For example, an inlet mesh with x, r, and 6 dimensions of 81 x 41 x 16, 
respectively, requires 53136 locations for the potential function alone. If a copy of the 
potential function at a previous step is required for extrapolation or other 
convergence acceleration p^rocedures, another 53136 locations are required. Any 
geometric information requires additional storage. Since these numbers are large 
compared to what is available, even for the vector computers such as the CYBER 203, 
it is essential that the scheme used for storing and accessing geometrical information 
be designed to minimize storage requirements. Disk storage may be used to provide 
more memory, but there are large time penalties for using disk storage; thus it is 
desirable to minimize its use. 

Storing a flag or index for each field point so that the code can determine if the field 
point is next to the surface, and if so, which surface point, requires at least one 
storage location for each field point. This type of scheme is conceptually simple, but 
inefficient. To minimize storage and maximize the percentage of vector code, an 
inverse scheme is used in this code. Surface-affected field points are processed by 
sweeping through the surface points. A bit array rather than a full word (64 bits) 
array is used to signal whether grid points are in the flowfield or not. A bit array is 
used for each of the three coordinates to mark whether field points are surface- 
adjacent in either or both directions in that coordinate. Total storage for the four bit 
arrays is 4/64 of 53136 or 3321 words. The bit arrays are used to inhibit storage of 
contributions to coefficients from regular field-point formulas when the field point is 
irregular (surface-affected). The correct contributions to coefficients of the 
finite-difference equations are then generated and stored by sweeping through the 
surface points. 

The geometric information stored for each surface point includes the index of the 
adjacent field point. 

Code efficiency on a vector machine is improved by ensuring that most computation is 
done in a vector mode. On the CYBER 203, longer vectors (up to a maximum length of 
65536) give higher computation rates. However, most vector computation involves 
intermediate results which are also vectors and must be stored. Hence there is a 
trade between the additional memory required for the intermediate vectors and the 
inceased computation rates obtained using longer vectors. For this problem, 
computing coefficients for one axial plane at one time appears to be a reasonable 
tradeoff. 

Surface points are grouped by type, adjacent in x, adjacent in r, adjacent in 0, or both 
a surface point and a field point, for each axial plane (plane perpendicular to axis). 
Quantities needed for coefficient calculation are gathered for all surface points of the 
same type, the coefficients are then calculated as a group using vector arithmetic, 
and the contributions to the coefficients are then scattered to the correct field-point 
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equations. This maximizes the number of vector calculations and minimizes the 
storage of geometric information. The primary penalty is the inefficiency of the 
gathering and scattering operations on the CYBER 203. The gather/scatter operations 
are more efficient on some other computers, and on the CYBER 205 which is the 
successor to the CYBER 203. 

2.7.2 PROCEDURE FOR SWEEPING 

The sequence of computations for taking a single relaxation sweep, or iteration, for a 
given mesh follows: 

I. The three components of velocity are calculated for all surface nodes. 

II. The field is swept one axial plane (constant x) at a time in the direction of 

increasing x. 

A. The three components of velocity at all field nodes are calculated for the 
axial plane ahead of the one being processed. (The velocities are saved for 
the plane being swept and the planes adjacent to the plane being swept.) 

B. The coefficients for the surface-node difference equations for all surface 
nodes on or adjacent to the plane being processed are calculated. 

C. The coefficients for the terms in the partial differential equation are 
calculated for all field nodes. 

D. The coefficients of the difference equation for all field-point nodes are 
calculated. 

E. The surface-node equations are used to eliminate the values of 0 at surface 
nodes from the field-point difference equations. 

F. The first radial line of difference equations is solved including the value for 
0 on the axis. 

G. The new value of 0 on the axis is used in the solution of all remaining 
radial lines. 

H. The new 0 values at field nodes are substituted into the surface-point 
difference equations and new values of 0 are obtained at all surface nodes. 

III. The entire field has been updated and a decision is made on whether to continue 

sweeping, extrapolate, change meshes, or stop. 

2.7.3 CONVERGENCE ACCELERATION 

The basic method described above is stable and convergent, but the rate of 
convergence for dense meshes can be very slow and computer times significant. There 
is substantial benefit from any procedure that will speed convergence and decrease 
the cost of a solution. Analysis and code development have been structured to 
facilitate the use of various techniques for convergence acceleration. One technique 
currently used is extrapolation of solutions as described in reference 6. The theory for 
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this technique can be obtained by approximating the solution process as a linear one 
and looking at the eigenvalues of a matrix. The number of eigenvalues is of the order 
of the number of field points. This approach behaves best when the eigenvalues are 
discrete, which is more likely when there are few of them. As a consequence, this 
technique works best for very coarse meshes and much less satisfactorily for fine 
meshes. 

If is <f) at iteration n, and II II is a norm, then 

X = ^ (51) 

■ ll/-^)|l 

and the extrapolation formula is 


9^. 


extrapolated 



(52) 


In practice, two norms are calculated after each sweep 



E 

all field 
nodes 





(53) 


and 



maximum 


k 


^i.j.k 


(54) 


where the maximum is taken over all field nodes. 

Extrapolation occurs whenever the eigenvalues and X 2 corresponding to II II ^ and 
II II 2 are steady and approximately the same for several sweeps. If they are 
reasonably close to being constant for several sweeps, it is assumed that the problem 
is behaving linearly and an extrapolation of the field will be of benefit. 

Another technique for convergence acceleration is initial convergence on a coarse 
mesh to obtain an approximate solution, continuation of the convergence on a 
medium-density mesh to obtain a more accurate approximation of the solution, and 
then final convergence on a fine mesh. The number of mesh for each coordinate is 
doubled for every successive mesh. There is approximately a factor of 16 change in 
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convergence rate per sweep, in terms of computer central-processor time, between 
successive meshes. This arises due to a factor of eight change in number of mesh 
points, and a factor of two in the rate of change per sweep. Once the solution nears 
convergence on one of the coarse meshes, there is no further gain in continuing 
relaxation on the mesh because of the large truncation errors. For the meshes used in 
the calculations of this report, there was about one-third decrease in run time using 
this approach. The reason for such a small improvement is the very slow convergence 
on the final mesh. It takes half as many sweeps to improve the solution from coarse 
meshes on the final mesh as it would to do the entire calculation using the finest 
mesh. Very good convergence on a coarse mesh takes less than 100 sweeps. Good' 
convergence using the fine mesh exclusively takes about 400 sweeps. Use of a still 
finer mesh starts to become very impractical due to the number of sweeps required 
and the time per sweep. 

Parametric cubic interpolation as described in appendix B is used to interpolate 
the solution from the current mesh to the next finer mesh. 

A further refinement on convergence acceleration is the multilevel technique 
described by Brandt^®\ This procedure uses coarse meshes to resolve the errors 
determined on finer meshes. The number of sweeps required on the finest mesh 
should be small (on the order of ten to thirty) and independent of the fineness of the 
mesh. The majority of sweeping would be on the coarsest mesh. As can be seen from 
the discussion of the preceding paragraph this technique becomes better and better, 
relative to the existing approach, as the meshes become denser. This analysis has 
been established with the anticipation of implementing a multilevel procedure in the 
future. 


2.8 GEOMETRY SPECIFICATION 

The basic analysis procedure described has no inherent limitations on kinds of 
geometries that can be analyzed, but any implementation on the computer does have 
limitations. Two sources of limitations are the number of mesh that can be used and 
the boundary conditions that can be applied. A limitation independent of the flow 
analysis code is the user's ability to specify the geometry and extract the 
surface/mesh intersections. 

The computer code for this analysis requires that all intersections of the mesh with 
the geometry, plus components of the surface normals, be input. The flow code then 
orders the points and makes some checks for completeness and consistency. The 
actual values could come from any source; no particular order is assumed by the flow 
code. The three-dimensional geometries used as samples in this report were described 
by a geometry system^®^ which breaks the surface into subsections or patches. The 
coordinates of the surface on the patch are specified by parametric bicubics. This 
leads to an explicit specification of the surface. A separate program exists to intersect 
the surface so described, with a computational mesh and to generate the coordinates 
and surface normals at the mesh intersections. 

2.9 RESULTS 

Comparisons between the analysis and experiment have been made for two 
asymmetric geometries. In addition, a mixer-lobe geometry has been analyzed to 
verify the program's general three-dimensional capabilities, although no 
experimental results are available for comparison. The first comparison between 
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theory and experiment is for an asymmetric tilt-nacelle V/STOL airplane inlet, tested 
in the NASA-Ames Research Center*s 40- by 80-ft wind tunnel. Figure 7 presents a 
graphic display of the inlet and spinner geometry. The geometry description used in 
the analysis is smooth and the entire inlet is specified. Only half the inlet is shown in 
figure 7, and curves are represented using linked straight-line segments due to 
limitations of the plotting equipment. The test and results of the test are described in 
reference 10. Tabular data from the test is found in reference 11. 

Predictions of the analysis have been compared with experimental measurements for 
three test cases (figs. 8, 9»and 10). There are two comparisons for an angle of attack 
of 60 degrees. The second comparison is for a greater airflow and freestream Mach 
number than the first. The last comparison is for an angle of attack of 90 degrees. 
The analysis predicts the experimentally measured Mach number distribution very 
closely. There is a slight, but consistent, under-prediction of the peak Mach number. 
This is possibly due to grid density and is discussed later. The only poor prediction is 
for the external lee side for 90 degree angle of attack (fig. 10). This probably is an 
interference or boundary-layer effect in the experiment. 



Figure 7. Graphical Display of V/STOL Airplane Inlet Geometry 
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Figure 8. Cowl Surface Mach Number Distribution for an Asymmetric VISTOL 
Airplane Inlet 
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Figure 9. Cowl Surface Mach Number Distribution for an Asymmetric VISTOL 
Airplane Inlet 
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Figure 10. Cowl Surface Mach Number Distribution for an Asymmetric V/STOL 
Airplane Inlet 





A series of comparisons have been made for a typical commercial-transport 
turbofan-engine inlet design. It has an asymmetric lip section varying both in section 
shape and contraction ratio. The contraction ratio on the lee side is 1.246 and on the 
windward side 1.28. The inlet was designed with an axisymmetric diffuser. The inlet 
centerline is tilted down five degrees with respect to the engine centerline and a 
short transition region exists between the inlet and the compressor (fan) face. This 
tilt of the inlet centerline is to reduce cruise drag by aligning the inlet with the local 
flowfield while still placing the engine and hence its thrust vector in the desired 
orientation. The inlet flowfield was computed using engine-centerline oriented 
coordinates, and as a consequence, the inlet appears very as 3 nnmetric. A cross-section 
of the inlet along with a computational mesh is shown in figure 11. Dimensions are in 
meters for the full-scale inlet. 

A 0.16 scale model of the inlet has been tested* in the Boeing 9- by 9-ft low-speed 
propulsion wind tunnel. The same inlet configuration in a 0.47 scale was tested^^^) 
the NASA- Ames Research Center’s 40- by 80-ft wind tunnel. The NAS A- Ames test 
model was modified slightly in the vicinity of the fan face to attach to a different 
engine than that for which the inlet was originally designed. The Boeing 9- by 9-ft 
tunnel test used a suction source instead of an engine at the fan face. The 9- by 9-ft 
tunnel test gave results at greater airflows than could be obtained in the Ames test. 

The comparisons between analysis and experiment for the 0.16 scale model are shown 
in figures 12 and 13. Figures 14 and 15 show results for the 0.47 scale model. The 
most significant discrepancy between the results of the analysis and the experimental 
measurements is for the greatest peak Mach number (fig. 13). Some of the 
discrepancy is due to blockage effects of a thick boundary layer in the diffuser on the 
windward side. This condition is near where the inlet boundary layer separates. The 
analysis curves are drawn to the compressor station. As the inlet is tilted relative to 
the engine, each curve ends at a different inlet centerline station. 

A single lobe of a mixer has been analyzed with this code. The lobe is shown in 
figure 16. This is a shaded-graphic representation^^^^ of the upper surface of the 
mixer. The inner surface was a constant-diameter cylinder. The lobe extended from 0 
degrees to 45 degrees and represents one lobe of an eight-lobe jet-engine mixer 
nozzle. Half a lobe could have been analyzed by taking advantage of the center plane 
of symmetry. The mixer shape was developed as a computer-program test case and 
not an actual mixer. It is shown to illustrate possible applications of the code. As 
there is no experiment with which to compare, the computed results are not shown. 
Use of the code for configurations other than inlets should include code validation by 
comparison of predictions with experiments for such configurations. 

2.10 ACCURACY 

The comparison between analysis and experiment for the V/STOL airplane inlet (figs. 
8, 9, and 10) show almost perfect agreement except at the point of peak Mach number 
on the windward side of the inlet. The progression of peak Mach number values for 
the case of figure 9 is 0.743, 0.893, and 1.007 for the three mesh levels used, clearly 
indicating a trend of an increase in peak Mach number for finer meshes. It is believed 
that the fine mesh used is the coarsest that will yield adequate accuracy for inlet 


*Boeing data, unpublished. 


35 


2.5 


2.0 

1.5 


1.0 


0.5 

ENGINE!^ 0 
INLETt 

0.5 


1.0 


1.5 

2.0 


2.5 



Figure 1 1. Computational Flowfield and Mesh in the Vicinity of a 
Commercial-Transport Turbofan-Engine Type Inlet 
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Figure 12. Cowl Surface Mach Number Distribution for an Asymmetric 
Commercial-Transport Turbofan-Engine Type Inlet 
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Figure 13. Cowl Surface Mach Number Distribution for an Asymmetric 
Commercial-Transport Turbofan-Engine Type Inlet 
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Figure 14. Cowl Surface Mach Number Distribution for an Asymmetric 
Commercial-Transport Turbofan-Engine Type Inlet 
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Figure 15. Cowl Surface Mach Number Distribution for an Asymmetric 
Commercial-Transport Turbofan-Engine Type Inlet 







design and analysis. The mesh used is an excellent tradeoff between cost and 
accuracy, since the solution for any significantly denser mesh would be extremely 
expensive to compute. It does appear that use of a finer mesh would result in higher 
predicted peak Mach numbers and slightly better agreement with experiment. Such 
calculations will be made when the multilevel technique of Brandt is implemented. 

There is an additional check on solution accuracy built into the program. The 
computed velocity and density profiles in the inlet are integrated at each axial 
station to compute mass flow, and the error relative to the mass flow enforced as a 
boundary condition is determined. Typical maximum errors are one-half to 
one-and-one-half percent. The magnitude of the error is a function of the complexity 
of the geometry, the mesh used, the degree of convergence, and the strength of any 
internal shocks. 

Overall, the program appears to predict inviscid flowfields to the accuracy that they 
can be measured experimentally. At flow conditions where Mach number gradients 
are large and boundary layers are thick, or shock-wave/boundary-layer interactions 
are present, agreement with experiment is not as good. It is expected that use of 
denser meshes and more stringent convergence tolerances would slightly improve 
agreement with experiment at the cost of significantly increased computer-time 
requirements. Improved accuracy for flows with significant viscous effects would 
require that viscous (boundary-layer) corrections be made to the inviscid flowfield. 

2.11 CONCLUSIONS 

The analysis technique described is capable of quite acceptable prediction for aircraft 
design and analysis purposes of flowfields for single-object geometries with 
essentially inviscid, irrotational flow. 

Reasonable predictions can be made for flows which are primarily irrotational except 
for thick, but attached, boundary layers. Geometry limitations are primarily due to 
the number of mesh lines the computer code can use and still have reasonable 
computation times. The restriction is that any major geometric variations have to 
have a scale of several mesh spacings. The restriction to inviscid irrotational flow is 
simply because viscous or boundary-layer effects have not been incorporated into the 
flow analysis. If these effects are significant, such as occurs with very thick or 
separated boundary layers, predictions will be poor. 
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3.0 USE OF THE PROGRAM 


This section discusses the use of the program exclusive of the generation of the 
geometry file. A detailed description of the input formats along with sample input 
and output files, and a troubleshooting guide, are presented. Also included is a 
discussion of selection of the computation mesh. 

The problem of specifying the geometry is not discussed. The specification of 
geometry and generation of the geometry part of the input file can be a major task. It 
has been treated as a separate topic with separate documentation. The code will 
accept geometry files from any source. The program does not require any special 
order to the points in the geometry file. 

The code, in its present form, can handle inlets with or without centerbodies, 
three-dimensional bodies, ducts, and bodies or inlets in a duct (e.g., wind tunnel). 
Specific limitations on geometries relative to mesh placement are discussed. 

3.1 COMPUTATIONAL PARAMETERS 

The analysis is programmed in CDC CYBER 200 FORTRAN language for the 

Control Data Corporation CYBER 203 computer. The code makes extensive use of the 
vector extensions of this language and is explicitly vectorized. Use is made of most of 
the million-word memory of the CYBER 203 at the NASA-Langley Research Center. 

A typical fine mesh for an inlet calculation is 69 axial-mesh by 37 radial-mesh by 16 
circumferential-mesh or 40848 points. The calculation of a solution, using just the 
fine mesh, to reasonable convergence takes approximately 400 sweeps and seven to 
eight minutes of central-processor time. Essentially the same solution can be 
obtained using a sequence of three grids in approximately five- and one-half minutes. 
A reasonable convergence pattern is 100 sweeps on the coarsest mesh (18 by 10 by 
16), 150 sweeps on the intermediate mesh (35 by 19 by 16), and 200 sweeps on the 
fine mesh. Typical times per sweep are 0.18 seconds per sweep on the coarse mesh, 
0.42 seconds per sweep on the intermediate mesh, and 1.1 seconds per sweep on the 
fine mesh. Some additional time is required for setting up geometry parameters, 
mesh changing, and printing the solution. Solutions have been calculated using a 
sequence of 4, 8, and 16 circumferential meshes, but the coarsest mesh had a AO of 
90 degrees and it is possible for this coarse spacing to cause some problems. 
Consequently, better convergence is presently obtained by keeping the number of 6 
mesh constant at 16. Run time is slightly dependent on the number of supersonic 
field points that have to be calculated. 

3-2 INPUT FORMAT 

The first two cards of the input deck are title cards and are printed at the start of the 
output for identification purposes. All input except the title cards is by means of 
order independent groups headed by keywords. The purposes of this particular input 
format are to allow certain groups to be optional, make the input file more readable, 
and to facilitate checking of input data by the program. Certain input groups are 
mandatory, and others are optional and may be omitted. All input, except for title 
cards, consists of numbers or words (depending on group) in fields of 10 columns wide, 
maximum of six fields per input line. All numbers are floating point and require a 
decimal point. Only the first four characters of keywords are checked. 
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There are certain interrelations among various input groups that have to be taken 
into account. If convergence is to be obtained on a sequence of meshes, the number of 
X, r, and 6 mesh can have only certain values. This is because coarser meshes are 
formed by deleting every other mesh line. Also, a compressor face, if there is one, 
must lie on an x mesh belonging to the coarsest mesh. There are restrictions on 
number of mesh and number of surface points relating to declared array lengths in 
the computer code. 

The program has the capability to use up to three mesh-density levels to provide 
more efficient convergence. The number of levels is controlled by the SWEEPS option. 
The mesh and geometry for the finest mesh level must be input. Coarser meshes for x 
and r are formed by deleting exactly every other mesh from the previous mesh. This 
places restrictions on the number of mesh allowed in the finest mesh, as the first and 
last mesh line have to remain when every other mesh is deleted. The 6 mesh is a 
special case. There is an option to control the manner in which the 0 mesh is varied 
between levels. The number of 0 mesh can be held constant for two successive levels, 
or every other d mesh value can be deleted for the coarser mesh. 

The program allows the use of planes of symmetry to cut the number of mesh needed 
to make a calculation. If the largest 6 mesh value input is 180.00 degrees, the plane 
0 degrees to 180 degrees is taken to be a plane of symmetry. If the largest 0 mesh 
input is less than 180 degrees the flow is assumed symmetrical about 0 degrees and 
the largest 0 value input. Zero degrees must always be input as a 0 mesh. 

3.2.1 INPUT-GROUP SUMMARY 


Keyword 

FREEstream 

XMESh 

RMESh 

TMESh 

GEOMetry 


Keyword 

COMPressor 

SWEEps 

THETa 


REQUIRED 

Description 

Speed of sound, freestream velocity, angle of attack, and 
angle of yaw 

Axial mesh values 

Radial mesh values 

Circumferential mesh values 

Surface/mesh intersections: coordinates and surface normals 
values 

OPTIONAL 

Description 

Indicates an inlet geometry and specifies inlet mass flow 
Convergence control parameters 

Control of number of d planes for each mesh-density level 
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Indicate special 0 differencing to be used 

Requests printout of various categories of geometrical 
information 

Requests surface flow variable printout at end of run 

Requests printout of flow variables along 0 constant cuts of 
surface 

Requests printout of flow variables at constant 6 cuts of 
flowfield 

Request for printout at other than level 3 for multilevel 
calculation 

OPTIONAL (Diagnostic) 

Description 

Requests print of coefficients, velocities and potential 
function for a specified axial cut and sweep number 

3.2.2 INPUT-GROUP DESCRIPTIONS 
FRj^STREAM 

This group specifies the velocity and orientation for the freestream relative to the 
geometry. 

The scaling of the velocities is essentially arbitrary except that they should be of 
order one to avoid difficulties with print formats. Note that 

Required input group, no default values. 


Card 1 

Cols. 

1-4 

'FREE’ 

Keyword 

Card 2 

Cols. 

1-10 

AINF 

aoo, freestream 
speed of sound 



11-20 

QINF 

q^o, freestream 
velocity 



21-30 

ALPHA 

a, angle of attack, 
degrees, arctan (yjua) 



31-40 

BETA 

P, angle of yaw. 


degrees, arctan 

Note: Input of 'TREE STREAM” which is 11 characters instead of "FREE” or 
"FREESTREAM” will draw an error message. 


SCDIff 
PRINt op 

SFLOw 
SURF pr 

FIELd pr 

IPRI 

Keyword 

DEBUg 
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XMES H 

RMESH 

TMESH 

These groups handle the input of the computational mesh, x, r, and theta, theta in 
degrees. The values do not have to be in any order. Theta mesh must include 0, 90, 180, 
and 270 degrees unless there is a plane of symmetry. Zero degrees must always be 
included. 

Required input groups, no default values. 

Card 1 Cols. 1-4 'XMES’ Keyword 

or 

'RMES’ 

or 

'TMES’ 

11-20 NX Number of mesh values to be read, 

or six per card. 

NR 

or 

NT 5«NX^ 101, 5^NR «81, 3=sNT:s41 

Card 2 Cols. 1-10 MESH(l) Axial, radial or circumferential 
11-20 MESH(2) location of mesh, 

21-30 six values per card, as many cards 

as required. Theta must 
be in degrees. 


Note: 

1) NX*NR*NT^56000 

2) NR*NT<800 

3) NX*NR^2800 

4) If three mesh levels are to be used in the analysis: 

NX-1=4,8,12,16 or 100. 

NR-1=4,8, 12,16 or 80 

5) The compressor face, if there is one, must lie on a x mesh in the coarsest mesh 
(i.e., XQp=Xi, i=l, or 5, or 9, or 13 etc. for a 3-level mesh). 

GEOM ETRY 

This group consists of the coordinates of the intersections of the mesh with the geometry 
and the direction cosines of the surface normal at each intersect. 
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Required input group, no default values. 


Card 1 

Cols. 

1-4 

'GEOM’ 

Keyword 




11-20 

NSURTOT 

number of intersects 

Card 2 

Cols. 

1-10 

SURFX(l) 

X 




11-20 

SURFR(l) 

r 

one 



21-30 

SURFT(l) 

6 (degrees) 

intersect 



31-40 

COSX(l) 

“x 

per 



41-50 

COSR(l) 

nr 

card 



51-60 

COST(l) 

rig 


Card 3 

Cols. 

1-10 

SURFX(2) 




Card 

NSURTOT+1 
COMP RESSOR (Optional 


Signals that there is an inlet geometry and specifies the Mach number at the compressor 
face. 

Card 1 Cols. 1-4 'COMP' Keyword 

11-20 AMACHCF Mach number to be enforced at 

compressor face. 

SWEE PS (Optional) 

This group controls the sweeping process by allowing control of the number of 
mesh-density levels, the maximum number of sweeps on each mesh-density level and a 
convergence criteria for each mesh-density level. 

Default is a three-level mesh with default values listed . 


Card 1 

Cols. 1-4 

'SWEE’ 

Ke3^ord 

Card 2 

Cols. 1-10 

NSWPM(l) 

Maximum number of sweeps 
on coarsest mesh (or total 
if one-level calculation) 


11-20 

NSWPM(2) 

Maximum number of sweeps 
on level 2, level 3 for a 
two-level calculation (zero 
implies a single-level calculation) 


21-30 

NSWPMC3) 

Maximum number of sweeps 
on level 3 (zero for a one- 
or two-level calculation) 
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Card 3 Cols. 1-10 CONVT(l) Convergence parameter, change in <f> 
11-10 CONVT(2) times 10®. Sweeping on a 

21-30 CONVT(3) level will stop when 

I X 10® ^ CONVT(i) 

CONVT(i) corresponds to the same level as NSWPM(i) 

Default values: 3 levels 

NSWPM(l) = 200. 

NSWPM(2) = 200. 

NSWPM(3) = 200. 

CONVT(l) = 1.0 
CONVT(2) = 1.0 
CONVT(3) = 1.0 

THETA (Optional) 

This group controls the number of $ grid used for each mesh-density level. The number 
of levels is controlled by either the SWEEPS option or the default value (three levels). 
The number of 6 intervals can be held constant or doubled for any mesh change. The 
values for the number of 0 grid must be consistent with the symmetry flag. 

Default for this option is no change in 9 grid for different levels. 


Card 1 

Cols. 

1-4 

•THET’ 

Ke3word 

Card 2 

Cols. 

1-10 

NTj 

Number of 0 grid for 
coarsest mesh-density level 



11-20 

NTg 

Number of 6 grid for level 2 
(level 3 for a two-level calculation) 



21-30 

NTg 

Number of 0 grid for level 3 
(for a three-level calculation) 


Note that if there is a plane of symmetry (less than 0 degrees to 360 degrees geometry 
input) 

NTi = NTj+i or NTj = (NTi+i+l)/2.0, 
otherwise, 

NTj = NTj+i or NTj = NTj+i/2.0 
SCDIff (Optional) 

This group allows use of special 9 differencing for improved accuracy with very coeirse 9 
meshes. 

Default is regular differencing. 
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Card 1 Cols. 


1-4 

11-20 


'SCDF 

ANUM 


Keyword 

0 . 0 Regular differencing, 

1.0 Special & differencing 
(any other value than 0.0 or 1,0 
will be treated as 0.0). 


PRINT OP (Optional) 

Inputs any or all of a group of ke 3 nvords to obtain printed output for certain geometrical 


quantities. 




Card 1 

Cols. 

1-4 TRIN’ 

Ke3nvord 

Card 2 

Cols. 

1-10 

Up to six ke 3 Twords 



11-20 

as described below. Can be 


in any order. 


SPIN PUT 

SPORD ER 

SPECPNTS 

TYPE 2 

MAP 


CUTS 


List of the surface points in the order read. 

List of surface points in the internal ordering used in the 
analysis. 

List of special points. 

List of Type 2 points. 

Lists of X, r, and B constant cuts of the surface. Lists 
include surface-point indexes, surface-point coordinates, 
arc length along the cuts, and components of the surface 
normals. 

Lists of X, r, and B constant cuts of the surface. Lists 
include surface-point indexes, surface-point coordinates 
and surface velocities. 


SFLQW (Optional) 

This group controls printing of flow properties along the surface. Default is printing of 
every fourth cut for all surfaces. This default corresponds to coarse-mesh cuts for a 
three-level calculation. 


Card 1 

Cols. 

1-4 

'SFLO’ 

Keyword 

Card 2 

Cols. 

1-10 

SKIPX 

See below. 



11-20 

SKIPR 




21-30 

SKIPT 



For the fine mesh, every constant x cut will be printed if SKIPX = 1. If SKIPX = 0 no 
cuts will be printed. Otherwise, cuts will be printed for x = X(I), I = 1,1 + SKIPX, 1 + 
2* SKIPX, etc. SKIPR and SKIIT' work the same for r and B constant cuts of the surfaces 
respectively. 



SURFACE PR (Optional) 


Surface properties along 6 constant cuts are printed for 0 equal to 0 degrees, 90 degrees, 
180 degrees, and 270 degrees for all runs except those with planes of symmetry. This 
option allows the printing of surface properties at other 0 mesh values. 


Card 1 

Cols. 

1-4 

•SURF’ 

Keyword 



11-20 

NSURPR 

Number of extra 6 printout 
planes. 

Card 2 

Cols. 

1-10 

THET(l) 

6 values for extra 



11-20 

THET(2) 

surface-property 
output. Six values 
per card (in degrees). 


FIEL D PR (Optional) 

This option determines for which 6 mesh values the field properties are to be printed. 


Card 1 Cols. 1-4 TIEL^ 
11-20 NFFPR 


Card 2 Cols. 1-10 THET(l) 
11-20 


Keyword 

Number of 6 values 
for which field 
properties are 
to be printed. 

Values of 6 mesh 
for printing field 
properties. Six values 
per card (in degrees). 


IPRI (Optional) 

This option allows printing of solution properties for the coarse meshes including 
mass-flow conservation computation. Default is no printout for coarse meshes. 

Card 1 Cols. 1-4 IPRP Ke 3 ^ord 

DEBUG (Optional) 

Diagnostic print option. Prints internal parameters for a given x mesh index and sweep 
number. 


Card 1 

Cols. 

1-4 

'DEBU’ 

Keyword 

Card 2 

Cols. 

1-10 

NPROPPR 

Sweep number 



11-20 

IPROPPR 

x-plane index 
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3.2.3 GEOMETRY SPECIFICATION 


Given a mesh, the user is required to provide all intersections between the mesh and the 
geometry to be analyzed. The information to be provided is the x, r, and d coordinates of 
the intersection, and the components of the unit normal to the surface at that 
intersection point, n^, nj.» and n^. The surface normals must be oriented such that the 
normals point from the surface into the flow. No particular order has to be followed in 
providing the intersections, but the set has to be complete. The program checks that the 
input points are complete and self-consistent, and generates internal maps connecting 
the points in x, r, and 6 constant cuts. 

The geometry points can come from any source. A package of procedures using 
parameteric bicubic patches has been developed to specify the geometry, and generate 
the geometry input file. It is documented separately (ref. 9). 

The program has been used to calculate the flow about three-dimensional ducts, inlets 
and bodies. Most of the program is extremely general and the flow about other 
geometries could possibly be calculated with minor program modifications. 

3.2.4 MESH SPECIFICATION 

The selection of the computational mesh can have a very significant effect on the cost 
and accuracy of a calculation. This is because the mesh affects the accuracy of the final 
converged solution through the truncation errors, and the rate of convergence of the 
solution. The influence of these effects is understood only in qualitative terms. In 
general, squarer meshes (Ax Ar === rAO) converge better. The critical mesh aspect ratio 
is Ax/rA0 since the use of line relsucation in the radial direction eliminates some of the 
effects of Ar. As the ratio of Ax/tAO departs significantly from unity, problems with 
convergence can occur. It is not possible, for inlet computations, to keep the far-field 
boundaries where they belong, keep a reasonable mesh near the inlet, and keep the mesh 
aspect ratio near unity everywhere. The best solution to this dilemma is to insure that 
any extreme aspect ratios occur in the far field where the solution changes little from the 
initial field. 

For accuracy, very fine meshes must be used in regions of the flow where large velocity 
gradients are expected. Inside ducts (including inlet ducts) it is desirable to have a 
reasonably dense mesh with mesh aspect ratios near unity in order to have good 
conservation of mass by the analysis. 

Table 3 gives a sample coarse mesh for a typical inlet configuration shown in figure 17. 
The fine meshes for x and r are generated by inserting three equally-spaced mesh lines 
between each of the coarse mesh lines. This will give a reasonable mesh for the inlet 
configuration shown. For vastly different geometries, the user may have to experiment 
with several possible meshes. It is very strongly suggested that the user make a sketch 
to scale, such as figure 17, showing the geometry and coarse mesh for any critical cuts, 
and take a good look at the sketch. 

Mesh Selection General Rules 

• Finer mesh in critical flow areas, that is areas where flow gradients are expected to 
be large. 
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xMESH 


r MESH 


6 MESH 


1. 

mSBM 

0.0 

2. 


0.3 R 

3. 


0.6 R 

4. 


0.8 R 

5. 


0.9 R 

6. 

■0.1 L 

1.0 R 

7. 

0.0 

1.2 R 

8. 

0.1 L 

1.6 R 


2.4 R 


10. 

0.4 L 

4.0 R 

11. 

0.6 L 


12. 

0.8 L 

10.0 R 

13. 

1.0 L 


14. 

1.2 L 


15. 

1.5 L 


16. 

2.0 L 



0.0 

22.5 

45.0 

67.0 

90.0 

112.5 

135.0 

157.5 

180.0 

202.5 

225.0 

247.5 

270.0 

292.5 

315.0 

337.5 








• Uniform or near-uniform meshes in ducts (including inlet ducts). 

• Mesh spacing should never be more than half or double between adjacent mesh (in 
the same coordinate). 

3.2.5 GEOMETRY LIMITATIONS 

This section is an attempt to give the user some understanding of what flowfields the 
code can successfully analyze, and those which it cannot. Also, it is intended to give the 
user some basic understanding of what the problems are when the program cannot 
successfully calculate the flowfield about a geometry as it may be possible to make a 
simple modification to the code to allow it to calculate a flowfield for which it previously 
was unsuccessful. 

The analysis is very general with respect to bodies lying internal to the computational 
volume. The primary limitations internal to the computational volume are mesh density 
and placement. Two difficulties may occur: the code may fail completely or it may give 
the wrong answer. Figure 18 illustrates limitations on mesh density in a passage 
between surfaces. If any one of a sequence of meshes is as shown in figure 18a the code 
will fail. A mesh such as that shown in figure 18b will not cause the code to fail, but 
accuracy probably would be poor. This may not be a problem if, for example, the mesh 
shown is the coarsest of a sequence of three meshes as the accuracy of the final solution 
is primarily dependent on the density of the finest mesh. Figure 19 illustrates that the 
mesh must be dense enough to resolve important features of the geometry. A solution 
calculated using the mesh shown in figure 19a would not predict any effects due to the 
wavy wall. A solution calculated using the mesh of figure 19b will predict at least some 
of the effects due to the wavy wall. 

The user has some control over placement of mesh, but certain mesh placement patterns 
could cause severe difficulties with convergence. Uniformly spaced mesh lines result in 
the best convergence behavior. A very irregular spacing of mesh lines, in an attempt to 
place the mesh where they are required for accuracy, may be self-defeating in that it 
might take a totally unreasonable number of sweeps to obtain convergence. This is a 
dilemma that would require numerical experiments to resolve for any given geometry. 

There are limitations on the total number of mesh that can be used which are built into 
the computer code. These could be changed by changing the code, but if the total number 
of mesh is increased, the run-time of the computer code also increases, and by a 
significantly greater factor. This run-time problem will be substantially alleviated when 
a version of the code with the multilevel procedure incorporated is available. 

Another limitation on the current version of the code is that there is no provision for 
application of a Kutta condition. While there is not believed to be any fundamental 
reason why a Kutta boundary-condition option could not be added to the code, it is 
thought to require a signficant engineering effort to incorporate. 

The external (edge of computational volume) boundary conditions allowed are those for 
an inlet, duct, body, or a body or inlet in a duct. Other geometries are possible with 
program modifications. The nature and degree of difficulty of such modifications would 
vary from problem to problem. 
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a.) INADEQUATE MESH FOR GEOMETRY RESOLUTION 



D.) ADEQUATE MESH FOR GEOMETRY RESOLUTION 


Figure 19. Mesh- Density Limitations Reiative to Surface Geometry 



The program uses cylindrical coordinates which are reasonably efficient for shapes that 
are nearly bodies of revolution. It is planned to incorporate a Cartesian-coordinate option 
in a later version of the code. Cartesian coordinates are much more efficient than 
cylindrical coordinates for analyzing some flowfields. 

3.2.6 SAMPLE INPUT FILES 

Several examples of input files (figs. 20a and 20b) are shown in an abbreviated form in 
order to help clarify the input file descriptions of Section 3.2.2. The files have been 
shortened by deleting most of the list of intersections and normals from the input files. 
The input groups on the files have been arranged for convenience of display by placing 
the "GEOMETRY” group as the last group on each file. 

The input file of figure 20a is the file used to generate the sample output of Section 3.3. 

3.3 PROGRAM OUTPUT 

The program output is arranged so that there is an introductory section, convergence 
history, and the flowfield solution. The introductory section (fig. 21a) presents program 
title, abstract, informative messages, and a list of consultants followed by specific 
information about the current run. The specific run information consists of run title, a 
repeat of program inputs such as free-stream conditions, convergence parameters, 
program options selected, etc. These are followed by the mesh. Input surface points and 
other geometry information can be printed at the user’s option. The headings used on 
surface point and geometry printouts are explained in Table 4. 

The convergence history (fig. 21b) includes sweep number, maximum and average 
residues for field points and surface points, location of maximum residue, number of 
supersonic points, the extrapolation parameters 1/(1-Xj^) and 1 /( 1 -X 2 )» flags to indicate 
convergence/divergence and extrapolation. Headings for the convergence history are 
listed in Table 5. 

Following the convergence history for each level, timing information for that level is 
printed (fig. 21c). 

After the timing information for the last level (and optionally for other levels), surface 
properties for 0 = 0 degrees, 90 degrees, 180 degrees and 270 degrees cuts (optionally for 
other 6 cuts), are printed (fig. 21d). The headings for surface-property printouts are listed 
in Table 6. Field properties may be printed at the user’s option (fig. 21e). Headings for 
field-property printouts are listed in Table 7. 

A mass-flow computation (fig. 2 If) is printed as a check on program convergence and/or 
accuracy. The mass flow for a duct is computed at each x station by integrating the 
computed solution profiles. The calculated mass flow is compared to the mass flow 
specified at the exit station. 

The last printout (fig. 21g) is optional (controlled by the SFLO input group) and allows 
printing of surface properties for x, r or ^ constant cuts. Headings are explained in 
Table 6. 
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XXX INLET APRIL 22,1980 IS IHETA PLANES 

175 KNOTS rtCF-O.SA ALPHA-25. 0 DEGREES R 13 C 3 

FREESTREAM 


1.0 

0.265 

25.0 

0.0 



SWEEPS 






100.0 

150.0 

200.0 

t 




J 

THETA LEV 

1 . 

1 • 




16.0 

16.0 

16.0 




CaMPRESSORO.64 





SFLOW 






4.0 

4.0 

4.0 




SCDIF 

1.0 





SURFACE PR2.0 





135.0 

225.0 





FIELD PR 

1.0 





180.0 






XMESH 

69.0 





-250.00000 

-218.75000 

-187.50000 

-156.25000 

-125.00000 

-108.75000 

-92.50000 

-76.25000 

-60.00000 

-50.00000 

-40.00000 

-30.00000 

-20.00000 

-12.50000 

-5.00000 

2.50000 

10.00000 

13.75000 

17.50000 

21.25000 

25.00000 

27.50000 

30.00000 

32.50000 

35.00000 

36.25000 

37.50000 

38,75000 

40.00000 

41.25000 

42.50000 

43.75000 

45.00000 

46.25000 

47.50000 

48.75000 

50.00000 

51.25000 

52.50000 

53.75000 

55.00000 

56.25000 

57.50000 

58.75000 

60.00000 

61.25000 

62.50000 

63.75000 

65.00000 

67.50000 

70.00000 

72.50000 

75.00000 

77.50000 

80.00000 

82.50000 

85.00000 

87.50000 

90.00000 

92.50000 

95.44000 

100.00000 

105.00000 

110.00000 

115.00000 

119.75000 

124.50000 

129.25000 

134.00000 




RMESH 

37.0 





0.00000 

5.00000 

10.00000 

15.00000 

20.00000 

23.75000 

27.50000 

31.25000 

35.00000 

37.50000 

40.00000 

42.50000 

45.00000 

48.25000 

47.50000 

48.75000 

50.00000 

52.50000 

55.00000 

57.50000 

60.00000 

62.50000 

65.00000 

67.50000 

70.00000 

77.50000 

85.00000 

92.50000 

100.00000 

112.50000 

125.00000 

137.50000 

150,00000 

175.00000 

200.00000 

225.00000 

250.00000 






TMESH 

16.0 





0.00000 

22.50000 

45.00000 

67.50000 

90.00000 

112.50000 

135.00000 

157.50000 

180.00000 

202.50000 

225.00000 

247.50000 

270.00000 

292.50000 

315.00000 

337.50000 



GEOMETRY 

1808. 





37.3789 

42.5000 

0.0000 

-.998896 

-.046984 

.000001 

38.7500 

44.6861 

22.5000 

-.825046 

.565030 

.006288 

38.7500 

44.6861 

337.5000 

-.825036 

.565045 

-.006288 

38.7500 

44.7495 

0.0000 

-.752486 

.658608 

.000000 

38.9909 

45.0000 

0.0000 

-.696756 

.717308 

-.000000 

38.9948 

45.0000 

22.5000 

-.745554 

.666445 

.000735 

38.9948 

45.0000 

337.5000 

-.745547 

.666453 

-.000735 

39,1819 

45,0000 

45.0000 

-.931933 

,361689 

,026095 

39.1819 

45.0000 

315.0000 

-.931928 

.361704 

-.026096 

40.0000 

45.0000 

64.4733 

-.996627 

-.026199 

.077768 

40.0000 

45.0000 

295.5275 

-.996627 

-.026183 

-.077769 

40.0000 

45.7338 

0.0000 

-.495294 

.868725 

-.000000 

40.0000 

45.8325 

22.5000 

-.530264 

.847771 

-.010236 

40.0000 

45.8325 

337.5000 

-.530260 

.847773 

.010235 

40.0000 

46.1618 

45.0000 

-.673944 

. 738626 

-.015237 

40.0000 

46.1618 

315.0000 

-.673939 

.738630 

.015236 

40.1000 

46.2500 

45.0000 

-.656499 

.754125 

-.017445 

40.1000 

46.2500 

315.0000 

-.656483 

.754139 

,017446 

40.3403 

46.2500 

67.5000 

-.924221 

.379650 

.041008 

40.3404 

46.2500 

292.5000 

-.924222 

.379648 

-.041009 

40.7701 

46.2500 

22.5000 

-.459877 

.887885 

-.013136 

40.7701 

46.2500 

337.5000 

-,459869 

.887889 

,013135 

41.0892 

46,2500 

0.0000 

-.375960 

.926636 

-.000000 


Figure 20a. Truncated Input File For An Inlet 
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MPY 23, 1980 


rORCED MIXER GEOMETRY 
M EXIT SET TO 0.45 
FREESTREPM 

1.0 0.45 0.0 

SHEEPS 

200.0 200.0 

1.0 1.0 

IPRI 

C0MPRESS0R0.45 
riELD PR 3.0 
0.0 7.5 15.0 

THETA LEV 

17.0 17.0 17.0 

XMESH 61.0 

0.0 0.1 0.2 

0.6 0.7 0.8 

1.2 1.3 1.4 

1.8 1.9 2.0 

2.4 2.5 2.6 

3.0 3.1 3.2 

3.6 3.7 3.8 

4.2 4,3 4.4 

4.8 4.9 5.0 

5.4 5.5 5.6 

6.0 

RMESH 37.0 

0.4 0.5 0.6 

1.0 1.1 1.2 

1.6 1.7 1.8 

2.2 2.3 2.4 

2.8 2.9 3.0 

3.4 3.5 3.6 

4.0 

TMESH 17.0 


0.0 0. 

9375 1 . 

875 

5.625 6. 

5625 7 . 

50 

11.25 12.1875 13 

GEOMETRY 3186. 

1.125 

0.0000 

2.9400 

0.0000 

0.0000 

2.9400 

.9375 

0.0000 

2.9400 

1 .8750 

0.0000 

2.9400 

2.8125 

0.0000 

2.9400 

3.7500 

0.0000 

2.9400 

4.6875 

0.0000 

2.9400 

5.6250 

0.0000 

2.9400 

6.5625 

0.0000 

2.9400 

7.5000 

0.0000 

2.9400 

8.4375 

0.0000 

2.9400 

9.3750 

0.0000 

2.9400 

10.3125 

0.0000 

2.9400 

11.2500 

0.0000 

2.9400 

12.1875 

0.0000 

2.9400 

13.1250 

0.0000 

2.9400 

14.0625 

0.0000 

2.9400 

15.0000 

.1000 

2.9396 

12.1875 

.1000 

2.9396 

13.1250 

,1000 

2.9396 

14.0625 

.1000 

2.9396 

15.0000 

.1000 

2.9397 

1 .8750 

.1000 

2.9397 

2.8125 

.1000 

2.9397 

3.7500 

.1000 

2.9397 

6.5625 

.1000 

2.9397 

9.3750 

.1000 

2.9397 

10.3125 


0.0 


0.3 

0.4 

0.5 

0.9 

1 .0 

1.1 

1.5 

1 .6 

1.7 

2.1 

2.2 

2.3 

2.7 

2.8 

2.9 

3.3 

3.4 

3.5 

3.9 

4.0 

4.1 

4.5 

4.6 

4.7 

5.1 

5.2 

5.3 

5.7 

5.8 

5.9 

0.7 

0.8 

0.9 

1.3 

1 .4 

1 .5 

1.9 

2.0 

2.1 

2.5 

2.6 

2.7 

3.1 

3.2 

3.3 

3.7 

3.8 

3.9 

2.8125 

3.75 

4.6875 

8.4375 

9.375 

10.3125 

14.0625 

15.0 



0.000000 

-1.000000 

0.000000 

.000025 

-1.000000 

0.000000 

.000031 

-1.000000 

0.000000 

-.000048 

-1.000000 

0.000000 

-.000018 

-1.000000 

0.000000 

-.000315 

-1.000000 

0.000000 

.000212 

-1.000000 

0.000000 

-.000398 

-1.000000 

0.000000 

.000271 

-1.000000 

0.000000 

.000197 

-1.000000 

0.000000 

-.000129 

-1.000000 

0.000000 

.000095 

-1.000000 

0.000000 

.000244 

-1.000000 

0.000000 

-.000085 

-1.000000 

0.000000 

-.000300 

-1 .000000 

0.000000 

-.000238 

-1.000000 

0.000000 

0.000000 

-1.000000 

0.000000 

-.007258 

-.999973 

.000800 

-.007751 

-.999969 

-.000977 

-.008033 

-.999968 

.000234 

-.008031 

-.999968 

0.000000 

-.005385 

-.999985 

-.000173 

-.005571 

-.999984 

-.000447 

-.005800 

-.999983 

.000246 

-.005002 

-.999987 

-.000605 

-.006119 

-.999981 

-.000975 

-.006129 

-.999979 

-.002120 


Figure 20b. Truncated Input File For A Mixer Lobe 
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P ^ 6 5 A - THRtE-ClH£NSi:NAL TRANSCNlC POTENTIAL FLOW PROGRAM 
VERSION OF MARCH 17* 1981 


RUN DATE - 03/23/81 


ABSTRACT - 


P465 COMPUTES SOLUTIONS FOR THE COMPLETE E^iUATlON FOR COMPRESSIBLE POTENTIAL FLOW. IT CAN PREDICT 
SUBSONIC OR TRANSONIC FLOy FIELDS ABOUT FULLY THREE-DIMENSIONAL INLETSt DUCTS AND BODIES* THE 
ANALYSIS USES CYLINDRICAL CCCROINATES* FINITE DIFFERENCES AND SUCCESSIVE LINE OVER RELAXATION CSLOR)* 
eXTRAPOLAUCN AND A SEQUENCE OF INCREASINGLY DENSER MESHES CAN BE USED TO ACCELERATE CONVERGENCE* 

THE ANALYSIS IS PROGRAMMED IK EXTENDED FORTRAN IV FOR THE CONTROL DATA CORPORATICN CYBER 203 
CCMPtTER. 

INITIAL RELEASE OF PA65 - VERSION A - NOVEMBER It 1980* 

REFERENCE - REYHNERt T* A** "TRANSONIC POTENTIAL FLCU COMPUTATION ABOUT THREE-DIMENSIONAL INLETSt 
DUCTS AND BOOlESt* AIAA PAPER 80-136A* SNOUHASSt COLORAOOt JULY 1980* 


PRCPRIETARY NOTICE - 


• THE COMPUTER PROGRAM* P465 - VERSION A* IS THE SOLE PROPERTY OF THE BOEING COMPANY UNTIL 

• NOVEMBER It 1983* DURING VHICH TIME NASA (THE NATIONAL AERONAUTICS AND SPACE ADHIN1STRAT1CN7 

• HAS RIGHTS OF USE. 


CONSULTATION - 


T* A* REYHNER 
R. G* JCRSTAO 
D* E. REUBUSH 


BOEING CCMHERCIAL AIRPLANE CO, 
BOEING COMPUTER SERVICES* INC. 
NASA LANGLEY RESEARCH CENTER 


(2061 237-2519 
(2061 656-5745 
(804> 827-2675 


RUN TITLE 


XXX INLET APRIL 22*1980 16 THETA PLANES 

175 KNOTS MCFrO,64 ALPHA=25*0 DEGREES R 19 C 3 


THE FLOW FIELD PARAMETERS ARE - 


AINF 

QINF 

ALPHA (ANGLE CF ATTACK) 
BETA (ANGLE OF YAH) 


l.OOOOQO 

*265000 

25.000 DEGREES 

•000 DEGREES 


6ECMETMY • 

INLET GEOMETRY - CCNPRESSCR FACE MACH NO* = .640 


MESH AND CONVERGENCE PARAMETERS - 


LEVEL NLMBER 


2 


3 


Figure 21a. Sample Output - Introductory Section 



NX 


69 


' 37 
16 
200 

1.0000 

SPECIAL SIN/COS DIFFERENCE QUOTIENTS USED FOR THETA COORDINATE 

SURFACE FLOU PROPERTIES PRINT REQUESTED - 

X CCASTANT CITS X=X(I) (I-l> OlUISIBLE BY 4 

R COKSTAAT CUTS R=R<JI (J-D DIVISIBLE BY 4 

THETA CCNSTAM CUTS THET A=THETAC K) <X*1> DIVISIBLE BY 4 


NR 

NT 

NAXIPUN NLPBEIi OF SHEEPS 
CONVERGENCE TEST VALUES «ao««6l 


18 

10 

16 

100 

1.0000 


35 

19 

16 

150 

1.0000 


Ih ADDITION rc THE 0, 90« 180« AND 270 DEGREE PLANES 
SURFACE POINT PRINTOUT HAS BEEN REQUESTED AT - 

135.0 DEG. 

225.0 CEE. 

PRINTOUT OF THE FLOyFIELO HAS BEEN REQUESTED AT - 

180.0 OEG. 


PESH - 


NX*NR*NT = 40848 


THIS DATA CASE USES 40848 CF THE AVAILABLE 56000 FIELD POINTS 


NR*NT 


592 the LIPIT IS 800 


X HESH 


1) 

•250.0000 

24) 

32.5000 

47) 

62.5000 

2) 

•218.7500 

25) 

35.0 00 0 

48) 

63.7500 

3) 

•187*5000 

261 

36.2500 

49) 

65.0000 

41 

•156.2500 

27) 

37.5000 

50) 

67.5Q00 

5) 

-125.0000 

281 

38.7500 

51) 

70.0QQQ 

6 ) 

•108.7500 

29) 

40.0000 

52) 

72.5000 

71 

-92.5000 

30) 

41.2500 

53) 

75.0000 

81 

-76.2500 

311 

42.5000 

54) 

77.5000 

91 

-60.0000 

32) 

43.7500 

55) 

80.0000 

10) 

-50.0000 

33) 

45.0000 

56) 

62.5000 

111 

-40.0000 

341 

46.2500 

57) 

85.0000 

121 

-30.0000 

35) 

47.5000 

58) 

87.5000 

131 

-20.0000 

36) 

48.7500 

59) 

90.0000 

14) 

-12.5000 

37) 

50.0000 

60) 

92.5000 

151 

-5.0000 

381 

51.2500 

611 

95.4400 

16) 

2.5000 

39) 

52.5000 

62) 

100 . 0000 

17) 

IC.OOOO 

40) 

53.7500 

63) 

105.0000 

181 

12.7500 

41 ) 

55.0000 

64 ) 

110.0000 

191 

17.5000 

42) 

56.2500 

65 ) 

115.0000 

20) 

?1 .?soo 

45 ) 

57.5000 

66) 

119.7500 


Figure 21a. Continued 



05 

21> 

2S .0000 

44) 

58. 7500 

67) 

124.5000 

to 

221 

27.5000 

45) 

60.0000 

68) 

129.2500 


23» 

30»QOGQ 

461 

61.2500 

69) 

134.0000 



« . * < 

r**.. R 

f*ESH 



1) 

.0000 

13) 

45.0000 

25) 

70.0000 

2) 

5.0000 

14) 

46.2500 

26) 

77.5000 

3) 

10.0000 

15) 

47.5000 

27) 

85.0000 

4) 

15.0000 

16) 

48.7500 

28) 

92.5000 

5) 

20.0000 

17) 

50.0000 

29) 

100.0000 

6) 

23.7500 

18) 

52.5000 

30) 

112.5000 

7) 

27.5000 

19) 

55.0000 

31) 

125.0000 

8) 

31.2500 

20) 

57.5000 

32) 

137-5000 

9) 

35.0000 

21 I 

60.0000 

33) 

150.0000 

10) 

37.5000 

22) 

62.5000 

34) 

175.0000 

11) 

40.0000 

231 

65.0000 

35) 

200.0000 

12) 

42.5000 

24) 

67.5000 

36) 

225.0000 





371 

250.0000 




• .* < 
*.* T 
. *« 1 

MESH 

t * * 


1) 

• 0000 

6) 

112.5000 

11) 

225.0000 

2) 

22.5000 

7) 

135.0000 

12) 

247.5000 

3) 

45.0000 

a> 

157 .5000 

13) 

270.0000 

4) 

67.5000 

9) 

180.0000 

14) 

292.5000 

5) 

90.0000 

10) 

202.5000 

IS) 

315.0000 





16) 

337.5000 

INPUT PROCESSING COHPLETED 






6E0METRT 

PROCESSING completed - 

LEVEL 3 



MESH FCR 

LEVEL 

2 - 

* * • < 

*• 

X 

* «* 

MESH 




1) 

•250.0000 

12) 

30.0000 

23) 

60.0000 


2) 

“187.5000 

13) 

35.0000 

24) 

62.5000 


3) 

-125.0000 

14) 

37.5C00 

25) 

65.0000 


4) 

-92.5000 

15) 

40.0000 

26) 

70.0000 


5) 

-60.0000 

16) 

42.5000 

27) 

75.0000 


6) 

-40.0000 

17) 

45.0C0O 

26) 

80.0000 


T) 

-20.0000 

18) 

47.5000 

29) 

85.0000 


8) 

-5.0000 

19) 

50.0000 

30) 

90.0000 


9) 

10.0000 

20) 

52.5000 

31) 

95.4400 


10) 

17.5000 

21) 

55.0000 

321 

105.0000 


11) 

25.0000 

22) 

57.5000 

33) 

115.0000 






34) 

124.5000 






35) 

134 .0000 



* * • i 

** 

R 

. * * 

MESH 

* 



Figure 21a. Continued 



05 

CO 


1) 

.0000 

7) 

45.0000 

13) 

70.0000 

2) 

10.0000 

8) 

47.5000 

14> 

85.0000 

3> 

20.0000 

9) 

50.0000 

15) 

100.0000 


27.5000 

10) 

55.0000 

16) 

125.0000 

5» 

35.0000 

11) 

60.0000 

17) 

150.0000 

6) 

40.0000 

12) 

65.0000 

lai 

200.0000 





19) 

250.0000 


T MESH 


1) 

.0000 

6) 

112.5000 

11) 

225.0000 

2) 

22.5000 

7) 

135.0000 

121 

247.5000 

3) 

45.0000 

8) 

157.5000 

13) 

270.0000 

4) 

67.5000 

9) 

180.0000 

14) 

292.5000 

5) 

90.0000 

10) 

202.5000 

15) 

315.0000 





16) 

337.5000 


SECMETPY PROCESSUG CCHPLETED - LEVEL 2 
MESH FCR LEVEL 1 - 


X MESH 


1) 

-250.0000 

7) 

35.0000 

13) 

65.0000 

2) 

-125.0000 

8) 

40.0000 

14) 

75.0000 

3) 

-60.0000 

9) 

45.0000 

15) 

85.0000 

4) 

-20.0000 

10) 

50.0000 

16) 

95.4400 

5) 

10.0000 

11) 

55.0000 

17) 

115.0000 

6) 

25.0000 

12) 

60.0000 

18) 

134.0000 


R MESH 


1) 

• 0000 

4) 

45.0000 

7) 

70.0000 

2) 

20.0000 

5) 

50.0000 

8) 

100.0000 

3) 

35.0000 

6) 

60.0000 

91 

150.0000 





10) 

250.0000 


T MESH 


1) 

• 0000 

6) 

112.5000 

11) 

225.0000 

2) 

22.5000 

7) 

135.QC00 

12) 

247.5000 

3) 

45.0000 

B) 

157.5000 

13) 

270.0000 

4) 

67.5000 

9) 

180.0000 

14) 

292.5000 

5) 

90.0000 

10) 

202.5000 

15) 

315.0000 





16) 

337.5000 


CECMETRY PROCESSING COMPLETED - LEVEL 1 


Figure 21a. Concluded 



CCilVERCENCE HISTORY 


LEVEL NUMBER 1 




.FiE-Ln pniMTs 



. Arr 

onruTC - 




SVEEP 

RVE RESICUE 

1 

J 

K 

MAX RESIDUE 

AVe RESIDUE 

I 

J 

K 

INDEX 

MAX RESIDUE 

M>1 

EIGENl 

1 

.77038E-fl3 

16 

4 

9 

.31551E-01 ** 

•51170E-O2 

16 

4 

9 

368 

•32044E-01 

0 

-.1650 7E-02 

2 

•E4244E-03 

15 

4 

9 

•16105E-01 • 

.29924E-02 

15 

4 

9 

320 

•16402E-01 

0 

•33528E«01 

3 

••9142E-03 

15 


9 

.11401E-01 * 

.24174E-02 

15 

4 

9 

32 0 

.11687E-01 

0 

.77278E401 

4 

•53283E-03 

15 

4 

9 

•89023E-02 * 

•19B64E-02 

15 

4 

9 

320 

.90558E-02 

0 

•75869E«01 

9 

•47503E-03 

16 

4 

9 

•71500E-02 * 

•16780E-02 

15 

4 

9 

320 

.723Q0e-02 

0 

•78805E«01 

6 

.43254E-03 

16 

4 

9 

•60413E-02 * 

•15232E-02 

15 

4 

9 

320 

•60862E-02 

0 

•11014E«02 

7 

«39265E*03 

15 

4 

9 

.50224E-02 * 

•13803E-02 

15 

4 

7 

314 

•50710E-02 

0 

•10744E«02 

8 

•36764E-03 

15 

4 

7 

•44443E-02 

•12890E-02 

15 

4 

7 

314 

•44841E-02 

0 

•15417E402 

9 

•33497E-03 

15 

4 

7 

.37364E-02 * 

•11630E-02 

15 

4 

7 

314 

«37691£-02 

0 

•10823C402 

10 

•31995E-03 

15 

4 

7 

•33970E-O2 

•10924E-02 

15 

4 

7 

314 

.34211E-02 

0 

•19646E402 

11 

•29861£*03 

15 

4 

7 

•29661E-02 * 

•10023E-02 

15 

4 

7 

314 

•29658E-02 

0 

•13768E4Q2 

12 

•28685E-03 

15 

4 

7 

•27457E-02 • 

•94837E-03 

15 

4 

7 

314 

•27605E-D2 

0 

•22353E402 

13 

•27118E*03 

15 

4 

7 

•24782E-02 * 

•d8278E-03 

15 

4 

7 

314 

•24908E-02 

0 

•16681E402 

14 

.26131E-C3 

15 

4 

7 

.23226E-02 • 

•64060E-03 

15 

4 

7 

314 

.23326E-02 

0 

•246B7E402 

15 

•24714E-03 

15 

4 

7 

.21069E-02 * 

•78332E-03 

15 

4 

7 

314 

•21158E-02 

0 

•16890E402 

16 

•23905E-03 

15 

4 

7 

.19938E-02 • 

.75024E-03 

15 

4 

7 

314 

•20009E-02 

0 

•27705E402 

17 

•23160E-03 

15 

4 

7 

•18930E-02 * 

•72134E-03 
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Figure 21b. Convergence History 
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Figure 21 d. Surface Properties - 6-Constant Cuts 
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Figure 21 e. Field Properties 
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Figure 21 f. Mass-Flow Conservation 
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Table 4. Headings for Surface Point and Surface Geometry Printout 


INDEX 

X 

R 

THETA 

NX 

NR 

NT 

S 

PHI,S-XC 
PHI,S-RC 
PHI,S-TC 
I ^ 

J ^ 

K , 

SURFARM 


TYPE 


Surface point index 

X 

r 

e 

n^, component of unit surface normal 
Uj., component of unit surface normal 
n^, component of unit surface normal 
s, arc length along cut 
(f>s^y velocity along x constant cut 
velocity along r constant cut 
velocity along 6 constant cut 

i, j, k of field point adjacent to surface 
point 

Ax, Ar or A9 between surface point and 
adjacent field point 

^ -1 = X intersect 

0 = theta intersect 

1 = r interset 

k 2 = both mesh and surface node 



Table 5. Convergence History Headings 


SWEEP 

**FIELD POINTS** 
AVE RESIDUE 


I 

"CONVERGING/ 

DIVERGING” 

** SURFACE POINTS** 


AVE RESIDUE 


INDEX I 

MAX RESIDUE / 


Relaxation sweep number, (n) 


Sum of I over all points 

i,j,k in the flowfield divided by 
and the number of such field points 

{ The indices i,j,k of the field point having 
the maximum change in 0 and the 
maximum value of 

divided by 

**, MAX RESIDUE decreasing, or 
***, MAX RESIDUE increasing 


Sum of over all surface points 

S divided by si^d the number of 

surface points 

{ The indices i,j,k of the field point next to 
the surface point with the maximum change 
in 0, the index of that surface point, 
and the maximum value of 1 1 
divided by 


M > 1.0 


number of field points for which Mach number >1.0 


EIGENl 


l/d-Xj), equations 51,53 


EIGEN2 


l/(l-\ 2 )> equations 51,54 


"EXTRAPOLATION FLAG” * indicates flowfield extrapolation using equation 

(52) was made after this sweep 
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Table 6. Surface-Point Printout Headings 


INDEX 

Surface point number 

X 

X 

R 

r 

THETA 

e 

S 

s, arc length along the cut of the surface 

MACH 

Mach Number 

CP 

Cp, coefficient of pressure, (p/poo-l)/^'yM2 
(if Moo = 0.0, then Cp = (p/p*-l)/Vfey, where 
the * indicates sonic conditions) 

PHI 


Q 

q = (u2 + Uj.2 + U^2)V<s 

PHI,S 

Ug, component of velocity along the cut 

U 

u 

U-RADIAL 


U-THETA 




Table 7. Field-Point Printout Headings 


R r 

MACH M 


PHI 

CP 


U 


0 

Cp, coefficient of pressure, (p/poo-l)/V 2 yM 2 
(if Moo = 0.0, then Cp = (p/p*-l)/y2'y, where the 
* indicates sonic conditions) 

u 


U-RADIAL 

U-THETA 


V 


V 


W 


w 


3.4 RUN QUALITY 

A major concern of the user is the degree to which the computed solution matches a 
real flowfield. This can he a relatively difficult question to answer. The primary 
effects are the nature of the flow analyzed, the accuracy of the analysis if a 
fully-converged solution is obtained, and the degree of convergence. The program 
predicts potential flow, which is inviscid irrotational flow. A necessary, but not 
sufficient, condition for acceptable flow prediction is that the flow to be analyzed not 
be greatly effected by viscous and rotational effects. 

If the flow is one for which potential flow is a good approximation, the next question 
is the acceptability of the mesh placement. If the flow is a duct, the mass-flow 
integration will give a check on accuracy. Otherwise, the best way to estimate 
accuracy (assuming a converged solution) is by comparison with experiments for the 
same or similar configurations. 

The question of convergence can be answered, to a large extent, by looking at the 
AVE RESIDUE in the convergence history. Values of 10‘® are usually satisfactory. 
For any given class of problems a few extra long runs can be made to obtain 10'^ or 
better convergence, and the results compared with lesser levels to see if there are any 
significant changes in the answers. 
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3.5 DIAGNOSTICS AND TROUBLESHOOTING 


Failure of the program to successfully provide a flow analysis for a data case could be 
due to a number of causes. Among these are input errors, errors in the program, and 
geometries the program was not written to handle. As obtaining help to determine if 
the program is in error can be very time consuming, it is preferable for a user to first 
check for input errors. The program does a considerable amount of input checking 
and many of the diagnostics are self-explanatory. The first step in searching for the 
cause of any run failure should be to check the printout for any diagnostic messages. 

One task of the program is checking the input geometry for completeness and 
self-consistency. Many of the error messages in this section are self-explanatory. One 
area that may be confusing is the messages from the subroutine which generates 
internal maps. There is a subroutine in the program that generates the maps that 
essentially list the surface points in sequence along x, r, and 0 constant cuts. If there 
are points missing, or points that are not consistent with each other, this subroutine 
will print diagnostic messages. The code works by starting at a surface point on a cut 
and then trying to follow the cut by finding the surface points in sequence on the cut. 

If the subroutine cannot find the next point on a cut, and the previous point was not 
a boundary point of the computational flowfield, an error message is printed and the 
points on the cut that have already been found will be printed. The problem area is at 
the end of the printed cut. Some more information may be obtained by looking at 
similar error printouts next to the one of interest as, usually, the parts of a cut on 
either side of a bad or missing surface point will be printed separately. The procedure 
for correcting the input error is to examine the area of difficulty, usually by sketches 
or plots, to determine the nature of the problem. If all surface points near the 
problem area are plotted versus the mesh, the problem is usually obvious. 

In interpreting the points printed out, the user should be aware that the mapping 
subroutine starts at a point and tracks the curve until it ends on a boundary or closes 
on itself. If the cut ends on a boundary, it then is followed in the opposite direction from 
the start point to the other end of the cut. If the subroutine is successful, the two 
segments of the curve are connected together in sequence. If the mapping routine fails 
while following the second segment, the points printed will be in sequence from the 
start point to the boundary point, followed by a sequence in the opposite sense from the 
point next to the start point to the problem area. 

The most common geometry problems are illustrated in figure 22. The problem of 
figure 22a is simply a point missing from the geometry file. The problem shown in 
figure 22b is inconsistent geometry points. Moving A to A' or B to B' would make the 
points consistent. This condition is usually caused by a tolerance problem in the 
geometry/mesh intersection program. Points A and B are very close to the mesh 
intersection and have been found by separate passes of the geometry processing code. 

A small error in position is enough to put the point on the wrong side of the mesh 
intersect. Often this problem can be fixed by replacing points A and B with a single 
point at the mesh intersect. Remember there is usually a point C in the third 
coordinate which should also be consolidated with A and B at the mesh intersect. 
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APPENDIX A 


TRIDIAGONAL EQUATION SOLVER 

The following is the logic for direct solution of a system of N simultaneous linear 
equations which are of tridiagonal form. If the unknown is <^j, j=l,2..,N, the 
equations are tridiagonal if they can be ordered such that the jth equations involves 
only 4 >- and 


The difference equations are: 





(A-1) 

Xj‘l'j-1 + Yj^j + Zj<^j+i = Gj 

2 ^ j ^ N-1 

(A-2) 

Xn'^N-i + Yn<^>n = G^. 


(A-3) 


For the above notation represents </>i j k as the lines for SLOR in this program are 
radial lines for which j is the radial index. 

The equations can be written in matrix notation as 

PW = R (A-4) 

where 


P = 


Yi Zi 

X2 Y2 Z2 

X3 Y3 Z3 


Xn-2 Yn-2 Zn-2 
Xn-1 Yn-1 
Xn 


Zn-1 

Yn 


(A-5) 


4>2 

^3 


G2 

G3 


w = 


(A-6) 


R = 


(A-7) 



Gn 
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mil II II I 


P can be decomposed into two triangular matrices, M and N, such that 
P= MN 
where 



Ml 

X2 M2 

Xg Mg 


Xn-2 Mn_2 

Xfj.l 

Xn Mn 


and 


N = 


Ni 

I N 2 
I Ng 


I Nn.2 

^ Nn.1 
I 


The following relations are then obtained; 


II 


Mj = Yj - XjNj.i 

2 j =sN 

Nj = Z/Mj 

1 j N-1 


(A-8) 


(A-9) 


(A-10) 


(A-11) 
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The matrix U is defined by 


MU = R 


where 


Ui 

U2 

U3 


u = 



(A-12) 


(A-13) 


which gives the relations 


Ui = Gi/Mi 

Uj = (Gj-XjUj.i)/Mj . 2 «j*N 


(A-14) 


When equations (A-8) and (A-12) are substituted into equation (A-4) and the resulting 
equation is multiplied by M'^, the following matrix equation is obtained : 


NW = U 


(A-15) 


which leads to the following equations for the cj)y 


- Un 


IsjsN.l 


(A-16) 


The solution procedure to obtain the <f>j is to use (A-11) to obtain Mj and Nj and (A-14) to 
obtain the Uj for all j. Then (A-16) is used to calculate the <t>y 
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APPENDIX B 


PARAMETRIC CUBIC INTERPOLATION 

Parametric cubic interpolation is used for interpolating to finer meshes as well as in 
several other areas of the program. It was chosen because odd-order interpolation 
leads to symmetric formulas, which are preferrable, and first-order or linear 
interpolation was thought to be too inaccurate. Also, cubic interpolation, as used with 
function and first derivatives, requires information at only two points. Given points A 
and B and a parameter u which is zero at A and 1 at B, then a function f(u) can be 
expressed as 

f(u) = Cq + CjU + C 2 U 2 + (B-1) 

however, the form 

f(u) = Fo(u)f(0) + Fi(u)f(l) + Do(u)fu(0) + Di(u)fjl) (B-2) 

where the Fj and are cubics in u, is preferrable for use in the code as it is easier to 
evaluate. 

It can be shown with a little algebra that 
Fo(u) = (1 -u)2(1 + 2u) 

Fi(u) = u2(3-2u) 

Dq(u) = u(l-u)2 

Di(u) - -u2(l-u). (B-3) 

The formulas are expressed in parametric form for simplicity. If x is the coordinate of 
interest and the known points are Xq and x^: 

u = (x-XqVCxi-Xq) (B-4) 

fu = (Xi-Xo)fx. (B-5) 


81 


I INI lllllllll III III III 



Uq = 0, Ui = 1.0, Vo = 0, Vi = 1.0, Wq = 0, Wi = 1.0. 

The above formula is used in the code for interpolation, except that the last four lines 
of equation B-6 are neglected. While this is incorrect, and introduces errors of 
unknown magnitude, it greatly simplifies the coding. Errors in the interpolation due 
to neglect of terms do not cause errors in the final flow solution, but instead increase 
the cost of a solution by requiring additional relaxation sweeps to eliminate the 
errors due to the interpolation process. 

Referring to figure B-1 for notation, the value of a function f at (u,v,w) can either be 
calculated from equation B-6 or by the following procedure: 




■t 


(0,0,1 

Figure B-1 Rectangular Box for Interpolation 
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First calculate, 


1 


fA = f(0,v,0) = 

S 1 

j=o ' 
1 

Fj(v)foj,o + Dj(v)fvo . J 


fuA = fu(0»V,0) = 

r I 

0=0 

1 



fwA = fw(0.v,0) = 

E 

j=0 

1 

[Fj(v)fwo j 0 + j_o] 

(B-7) 

^wua ~ ~ 

E 

j=o 

[Fj(v)fwuoj_o + ^j('')^uvwoj,o] 



The same quantities are calculated at A', A" and A'" by similar formulas. Then 
calculate 

fg = f(u,v,0) = Fo(u)fA + Fi(u)fA^ + Do(u)fuA + Di(u)f^^, 

fwB = fw(ii»v,0) = Fo(u)f^^ + Fj(u)f^^, + Do(u)f^^ + Di(u)f^u^,. (B-8) 

fgr and are calculated similarly 

Then 

fc = f(u,v,w) = Fo(w)fg + Fi(w)fg. +Do(w)f^^ + Di(w)f^g,. (B-9) 

The above procedure is exactly equivalent to using equation B-6. The advantage of 
the above procedure is that it can easily be used with minor modifications if one or 
more corners of the box are missing, which can occur in the vicinity of surfaces. This 
gives a workable interpolation procedure near surfaces, however, the results, when a 
corner of the box is missing, are not unique. That is, interpolating first in v, then in 
u, and finally in w as described, is not equivalent to another sequence. It is 
equivalent only if the rectangular box is complete. 
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